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: example/polynomial.py
===================================================================
--- example/polynomial.py	(revision d19962c0471a0eee538ce03020b5b226ba2752f6)
+++ example/polynomial.py	(revision d19962c0471a0eee538ce03020b5b226ba2752f6)
@@ -0,0 +1,12 @@
+from numpy import inf
+parameters = [
+    ["n", "", 1, [1,5], "", "number of coefficients (or degree+1)"],
+    ["c[n]", "", 0, [-inf, inf], "", "coefficients to c_n x^n"],
+]
+
+Iq = r"""
+    int int_n = (int)n;
+    double result = c[int_n-1];
+    for (int k=int_n-2; k >= 0; k--) { result = result*q + c[k]; }
+    return result;
+    """
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
+++ sasmodels/bumps_model.py	(revision 21b116f73d822594145622e342e0308fc3ecd14a)
@@ -81,23 +81,20 @@
     from bumps.names import Parameter
 
-    pars = {}
-    for p in model_info['parameters']:
+    pars = {}     # floating point parameters
+    pd_types = {} # distribution names
+    for p in model_info['parameters'].call_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 xpart, xdefault, xlimits in [
-                ('_pd', 0., pars[name].limits),
-                ('_pd_n', 35., (0, 1000)),
-                ('_pd_nsigma', 3., (0, 10)),
-            ]:
-            xname = name + xpart
-            xvalue = kwargs.pop(xname, xdefault)
-            pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits)
-
-    pd_types = {}
-    for name in model_info['partype']['pd-2d']:
-        xname = name + '_pd_type'
-        xvalue = kwargs.pop(xname, 'gaussian')
-        pd_types[xname] = xvalue
+        if p.polydisperse:
+            for part, default, limits in [
+                    ('_pd', 0., pars[p.name].limits),
+                    ('_pd_n', 35., (0, 1000)),
+                    ('_pd_nsigma', 3., (0, 10)),
+                ]:
+                name = p.name + part
+                value = kwargs.pop(name, default)
+                pars[name] = Parameter.default(value, name=name, limits=limits)
+            name = p.name + '_pd_type'
+            pd_types[name] = kwargs.pop(name, 'gaussian')
 
     if kwargs:  # args not corresponding to parameters
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision f2473145a6baa12aa3ab15421a80d2f5d004935c)
+++ sasmodels/compare.py	(revision 8bd7b77da137d537b08a21d51f0a8000689d4cc6)
@@ -38,5 +38,4 @@
 from . import core
 from . import kerneldll
-from . import product
 from .data import plot_theory, empty_data1D, empty_data2D
 from .direct_model import DirectModel
@@ -210,4 +209,5 @@
     Choose a parameter range based on parameter name and initial value.
     """
+    # process the polydispersity options
     if p.endswith('_pd_n'):
         return [0, 100]
@@ -222,17 +222,13 @@
         else:
             return [-180, 180]
+    elif p.endswith('_pd'):
+        return [0, 1]
     elif 'sld' in p:
         return [-0.5, 10]
-    elif p.endswith('_pd'):
-        return [0, 1]
     elif p == 'background':
         return [0, 10]
     elif p == 'scale':
         return [0, 1e3]
-    elif p == 'case_num':
-        # RPA hack
-        return [0, 10]
     elif v < 0:
-        # Kxy parameters in rpa model can be negative
         return [2*v, -2*v]
     else:
@@ -240,5 +236,5 @@
 
 
-def _randomize_one(p, v):
+def _randomize_one(model_info, p, v):
     """
     Randomize a single parameter.
@@ -246,9 +242,24 @@
     if any(p.endswith(s) for s in ('_pd_n', '_pd_nsigma', '_pd_type')):
         return v
+
+    # Find the parameter definition
+    for par in model_info['parameters'].call_parameters:
+        if par.name == p:
+            break
     else:
-        return np.random.uniform(*parameter_range(p, v))
-
-
-def randomize_pars(pars, seed=None):
+        raise ValueError("unknown parameter %r"%p)
+
+    if len(par.limits) > 2:  # choice list
+        return np.random.randint(len(par.limits))
+
+    limits = par.limits
+    if not np.isfinite(limits).all():
+        low, high = parameter_range(p, v)
+        limits = (max(limits[0], low), min(limits[1], high))
+
+    return np.random.uniform(*limits)
+
+
+def randomize_pars(model_info, pars, seed=None):
     """
     Generate random values for all of the parameters.
@@ -261,5 +272,5 @@
     with push_seed(seed):
         # Note: the sort guarantees order `of calls to random number generator
-        pars = dict((p, _randomize_one(p, v))
+        pars = dict((p, _randomize_one(model_info, p, v))
                     for p, v in sorted(pars.items()))
     return pars
@@ -294,10 +305,10 @@
         # Make sure phi sums to 1.0
         if pars['case_num'] < 2:
-            pars['Phia'] = 0.
-            pars['Phib'] = 0.
+            pars['Phi1'] = 0.
+            pars['Phi2'] = 0.
         elif pars['case_num'] < 5:
-            pars['Phia'] = 0.
-        total = sum(pars['Phi'+c] for c in 'abcd')
-        for c in 'abcd':
+            pars['Phi1'] = 0.
+        total = sum(pars['Phi'+c] for c in '1234')
+        for c in '1234':
             pars['Phi'+c] /= total
 
@@ -306,19 +317,13 @@
     Format the parameter list for printing.
     """
-    if is2d:
-        exclude = lambda n: False
-    else:
-        partype = model_info['partype']
-        par1d = set(partype['fixed-1d']+partype['pd-1d'])
-        exclude = lambda n: n not in par1d
     lines = []
-    for p in model_info['parameters']:
-        if exclude(p.name): continue
+    parameters = model_info['parameters']
+    for p in parameters.user_parameters(pars, is2d):
         fields = dict(
-            value=pars.get(p.name, p.default),
-            pd=pars.get(p.name+"_pd", 0.),
-            n=int(pars.get(p.name+"_pd_n", 0)),
-            nsigma=pars.get(p.name+"_pd_nsgima", 3.),
-            type=pars.get(p.name+"_pd_type", 'gaussian'))
+            value=pars.get(p.id, p.default),
+            pd=pars.get(p.id+"_pd", 0.),
+            n=int(pars.get(p.id+"_pd_n", 0)),
+            nsigma=pars.get(p.id+"_pd_nsgima", 3.),
+            type=pars.get(p.id+"_pd_type", 'gaussian'))
         lines.append(_format_par(p.name, **fields))
     return "\n".join(lines)
@@ -454,5 +459,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
@@ -661,13 +667,17 @@
     """
     # Get the default values for the parameters
-    pars = dict((p.name, p.default) for p in model_info['parameters'])
-
-    # Fill in default values for the polydispersity parameters
-    for p in model_info['parameters']:
-        if p.type in ('volume', 'orientation'):
-            pars[p.name+'_pd'] = 0.0
-            pars[p.name+'_pd_n'] = 0
-            pars[p.name+'_pd_nsigma'] = 3.0
-            pars[p.name+'_pd_type'] = "gaussian"
+    pars = {}
+    for p in model_info['parameters'].call_parameters:
+        parts = [('', p.default)]
+        if p.polydisperse:
+            parts.append(('_pd', 0.0))
+            parts.append(('_pd_n', 0))
+            parts.append(('_pd_nsigma', 3.0))
+            parts.append(('_pd_type', "gaussian"))
+        for ext, val in parts:
+            if p.length > 1:
+                dict(("%s%d%s"%(p.id,k,ext), val) for k in range(p.length))
+            else:
+                pars[p.id+ext] = val
 
     # Plug in values given in demo
@@ -774,8 +784,8 @@
 
     if len(engines) == 0:
-        engines.extend(['single', 'sasview'])
+        engines.extend(['single', 'double'])
     elif len(engines) == 1:
-        if engines[0][0] != 'sasview':
-            engines.append('sasview')
+        if engines[0][0] != 'double':
+            engines.append('double')
         else:
             engines.append('single')
@@ -807,5 +817,5 @@
     #pars.update(set_pars)  # set value before random to control range
     if opts['seed'] > -1:
-        pars = randomize_pars(pars, seed=opts['seed'])
+        pars = randomize_pars(model_info, pars, seed=opts['seed'])
         print("Randomize using -random=%i"%opts['seed'])
     if opts['mono']:
@@ -879,12 +889,12 @@
         model_info = opts['def']
         pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars'])
+        # Initialize parameter ranges, fixing the 2D parameters for 1D data.
         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'].user_parameters(is2d=False):
+                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/compare_many.py
===================================================================
--- sasmodels/compare_many.py	(revision ce346b6e747be779612d9528b60d52c717c995af)
+++ sasmodels/compare_many.py	(revision 8bd7b77da137d537b08a21d51f0a8000689d4cc6)
@@ -20,5 +20,4 @@
 
 from . import core
-from . import generate
 from .compare import (MODELS, randomize_pars, suppress_pd, make_data,
                       make_engine, get_pars, columnize,
@@ -109,5 +108,5 @@
 
     if is_2d:
-        if not model_info['has_2d']:
+        if not model_info['parameters'].has_2d:
             print(',"1-D only"')
             return
@@ -125,7 +124,5 @@
         try:
             result = fn(**pars)
-        except KeyboardInterrupt:
-            raise
-        except:
+        except Exception:
             traceback.print_exc()
             print("when comparing %s for %d"%(name, seed))
@@ -246,5 +243,5 @@
         base = sys.argv[5]
         comp = sys.argv[6] if len(sys.argv) > 6 else "sasview"
-    except:
+    except Exception:
         traceback.print_exc()
         print_usage()
Index: sasmodels/convert.json
===================================================================
--- sasmodels/convert.json	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/convert.json	(revision 6e7ff6d026d856e1eb236a05975a2600ba9ba2c9)
@@ -602,4 +602,9 @@
     "RPAModel", 
     {
+      "N1": "Na", "N2": "Nb", "N3": "Nc", "N4": "Nd",
+      "L1": "La", "L2": "Lb", "L3": "Lc", "L4": "Ld",
+      "v1": "va", "v2": "vb", "v3": "vc", "v4": "vd",
+      "b1": "ba", "b2": "bb", "b3": "bc", "b4": "bd",
+      "Phi1": "Phia", "Phi2": "Phib", "Phi3": "Phic", "Phi4": "Phid",
       "case_num": "lcase_n"
     }
@@ -620,7 +625,16 @@
     }
   ], 
+  "_spherepy": [
+    "SphereModel",
+    {
+      "sld": "sldSph",
+      "radius": "radius",
+      "sld_solvent": "sldSolv"
+    }
+  ],
   "spherical_sld": [
     "SphereSLDModel", 
     {
+      "n": "n_shells",
       "radius_core": "rad_core", 
       "sld_solvent": "sld_solv"
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision f2473145a6baa12aa3ab15421a80d2f5d004935c)
+++ sasmodels/convert.py	(revision 8bd7b77da137d537b08a21d51f0a8000689d4cc6)
@@ -206,5 +206,5 @@
         elif name == 'rpa':
             # convert scattering lengths from femtometers to centimeters
-            for p in "La", "Lb", "Lc", "Ld":
+            for p in "L1", "L2", "L3", "L4":
                 if p in oldpars: oldpars[p] *= 1e-13
         elif name == 'core_shell_parallelepiped':
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 4d767110d389aa65268b4b929ac8df1413460380)
+++ sasmodels/core.py	(revision ee8f73486ed89c512902019a13b040e576091f49)
@@ -24,5 +24,5 @@
     from . import kernelcl
     HAVE_OPENCL = True
-except:
+except Exception:
     HAVE_OPENCL = False
 
@@ -64,5 +64,5 @@
     """
     try: s + ''
-    except: return False
+    except Exception: return False
     return True
 
@@ -170,5 +170,5 @@
 
 
-def get_weights(model_info, pars, name):
+def get_weights(parameter, values):
     """
     Generate the distribution for parameter *name* given the parameter values
@@ -178,11 +178,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 = float(values.get(parameter.name, parameter.default))
+    relative = parameter.relative_pd
+    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)
+    if npts == 0 or width == 0:
+        return [value], []
     value, weight = weights.get_weights(
         disperser, npts, width, nsigma, value, limits, relative)
@@ -198,4 +200,5 @@
     """
     value, weight = zip(*pars)
+    weight = [w if w else [1.] for w in weight]
     value = [v.flatten() for v in meshgrid(*value)]
     weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
@@ -216,12 +219,32 @@
     *mono* is True if polydispersity should be set to none on all parameters.
     """
-    fixed_pars = [pars.get(name, kernel.info['defaults'][name])
-                  for name in kernel.fixed_pars]
+    parameters = kernel.info['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':
+        active = lambda name: name in parameters.pd_1d
+    elif kernel.dim == '2d':
+        active = lambda name: name in parameters.pd_2d
+    else:
+        active = lambda name: True
+
+    vw_pairs = [(get_weights(p, pars) if active(p.name)
+                 else ([pars.get(p.name, p.default)], []))
+                for p in parameters.call_parameters]
+
+    details, weights, values = build_details(kernel, vw_pairs)
+    return kernel(details, weights, values, cutoff)
+
+def build_details(kernel, pairs):
+    values, weights = zip(*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 details, weights, values
+
 
 def call_ER_VR(model_info, vol_pars):
@@ -256,9 +279,9 @@
         return 1.0
     else:
-        vol_pars = [get_weights(model_info, values, name)
-                    for name in model_info['partype']['volume']]
+        vol_pars = [get_weights(parameter, values)
+                    for parameter in model_info['parameters'].call_parameters
+                    if parameter.type == 'volume']
         value, weight = dispersion_mesh(vol_pars)
         individual_radii = ER(*value)
-        #print(values[0].shape, weights.shape, fv.shape)
         return np.sum(weight*individual_radii) / np.sum(weight)
 
@@ -273,6 +296,7 @@
         return 1.0
     else:
-        vol_pars = [get_weights(model_info, values, name)
-                    for name in model_info['partype']['volume']]
+        vol_pars = [get_weights(parameter, values)
+                    for parameter in model_info['parameters'].call_parameters
+                    if parameter.type == 'volume']
         value, weight = dispersion_mesh(vol_pars)
         whole, part = VR(*value)
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision d6f5da6a131a8de1b3e3d9a12af714a8d32aa4be)
+++ sasmodels/data.py	(revision ee8f73486ed89c512902019a13b040e576091f49)
@@ -358,7 +358,5 @@
         try:
             return fn(*args, **kw)
-        except KeyboardInterrupt:
-            raise
-        except:
+        except Exception:
             traceback.print_exc()
 
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision 4d767110d389aa65268b4b929ac8df1413460380)
+++ sasmodels/direct_model.py	(revision ee8f73486ed89c512902019a13b040e576091f49)
@@ -67,6 +67,4 @@
             self.data_type = 'Iq'
 
-        partype = model.info['partype']
-
         if self.data_type == 'sesans':
             
@@ -82,6 +80,6 @@
             q_mono = sesans.make_all_q(data)
         elif self.data_type == 'Iqxy':
-            if not partype['orientation'] and not partype['magnetic']:
-                raise ValueError("not 2D without orientation or magnetic parameters")
+            #if not model.info['parameters'].has_2d:
+            #    raise ValueError("not 2D without orientation or magnetic parameters")
             q = np.sqrt(data.qx_data**2 + data.qy_data**2)
             qmin = getattr(data, 'qmin', 1e-16)
@@ -172,5 +170,5 @@
     def _calc_theory(self, pars, cutoff=0.0):
         if self._kernel is None:
-            self._kernel = self._model.make_kernel(self._kernel_inputs)  # pylint: disable=attribute-dedata_type
+            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)
@@ -243,5 +241,5 @@
         try:
             values = [float(v) for v in call.split(',')]
-        except:
+        except Exception:
             values = []
         if len(values) == 1:
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
+++ sasmodels/generate.py	(revision 6e7ff6d026d856e1eb236a05975a2600ba9ba2c9)
@@ -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:
 
@@ -156,25 +125,4 @@
 An *model_info* dictionary is constructed from the kernel meta data and
 returned to the caller.
-
-The model evaluator, function call sequence consists of q inputs and the return vector,
-followed by the loop value/weight vector, followed by the values for
-the non-polydisperse parameters, followed by the lengths of the
-polydispersity loops.  To construct the call for 1D models, the
-categories *fixed-1d* and *pd-1d* list the names of the parameters
-of the non-polydisperse and the polydisperse parameters respectively.
-Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models.
-The *pd-rel* category is a set of those parameters which give
-polydispersitiy as a portion of the value (so a 10% length dispersity
-would use a polydispersity value of 0.1) rather than absolute
-dispersity such as an angle plus or minus 15 degrees.
-
-The *volume* category lists the volume parameters in order for calls
-to volume within the kernel (used for volume normalization) and for
-calls to ER and VR for effective radius and volume ratio respectively.
-
-The *orientation* and *magnetic* categories list the orientation and
-magnetic parameters.  These are used by the sasview interface.  The
-blank category is for parameters such as scale which don't have any
-other marking.
 
 The doc string at the start of the kernel module will be used to
@@ -204,24 +152,19 @@
 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
 
+from .modelinfo import ModelInfo, Parameter, make_parameter_table, set_demo
 from .custom import load_custom_kernel_module
 
-PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description']
-Parameter = namedtuple('Parameter', PARAMETER_FIELDS)
-
-C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
+TEMPLATE_ROOT = dirname(__file__)
 
 F16 = np.dtype('float16')
@@ -232,10 +175,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
@@ -331,4 +268,5 @@
     raise ValueError("%r not found in %s" % (filename, search_path))
 
+
 def model_sources(model_info):
     """
@@ -339,17 +277,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):
@@ -362,5 +297,5 @@
     if dtype == F16:
         fbytes = 2
-        source = _F16_PRAGMA + _convert_type(source, "half", "f")
+        source = _convert_type(source, "half", "f")
     elif dtype == F32:
         fbytes = 4
@@ -368,5 +303,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
@@ -411,28 +346,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_function_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):
     """
@@ -453,197 +432,246 @@
     # 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
+
+    # 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 = partable.form_volume_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] + partable.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,
-        }
-
-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
-
-def process_parameters(model_info):
-    """
-    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)
-    if model_info.get('demo', None) is None:
-        model_info['demo'] = model_info['defaults']
-    model_info['has_2d'] = partype['orientation'] or partype['magnetic']
-
+        pars = [qx, qy] + partable.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 partable.kernel_parameters))
+
+    # Define the function calls
+    if partable.form_volume_parameters:
+        refs = _call_pars("_v.", partable.form_volume_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) 1.0"
+    source.append(call_volume)
+
+    refs = ["_q[_i]"] + _call_pars("_v.", partable.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.", partable.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"%partable.max_pd)
+    source.append("#define NPARS %d"%partable.npars)
+
+    # 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)
+
+class CoordinationDetails(object):
+    def __init__(self, model_info):
+        parameters = model_info['parameters']
+        max_pd = parameters.max_pd
+        npars = parameters.npars
+        par_offset = 4*max_pd
+        self.details = np.zeros(par_offset + 3*npars + 4, 'i4')
+
+        # generate views on different parts of the array
+        self._pd_par     = self.details[0*max_pd:1*max_pd]
+        self._pd_length  = self.details[1*max_pd:2*max_pd]
+        self._pd_offset  = self.details[2*max_pd:3*max_pd]
+        self._pd_stride  = self.details[3*max_pd:4*max_pd]
+        self._par_offset = self.details[par_offset+0*npars:par_offset+1*npars]
+        self._par_coord  = self.details[par_offset+1*npars:par_offset+2*npars]
+        self._pd_coord   = self.details[par_offset+2*npars:par_offset+3*npars]
+
+        # theta_par is fixed
+        self.details[-1] = parameters.theta_offset
+
+    @property
+    def ctypes(self): return self.details.ctypes
+
+    @property
+    def pd_par(self): return self._pd_par
+
+    @property
+    def pd_length(self): return self._pd_length
+
+    @property
+    def pd_offset(self): return self._pd_offset
+
+    @property
+    def pd_stride(self): return self._pd_stride
+
+    @property
+    def pd_coord(self): return self._pd_coord
+
+    @property
+    def par_coord(self): return self._par_coord
+
+    @property
+    def par_offset(self): return self._par_offset
+
+    @property
+    def num_active(self): return self.details[-4]
+    @num_active.setter
+    def num_active(self, v): self.details[-4] = v
+
+    @property
+    def total_pd(self): return self.details[-3]
+    @total_pd.setter
+    def total_pd(self, v): self.details[-3] = v
+
+    @property
+    def num_coord(self): return self.details[-2]
+    @num_coord.setter
+    def num_coord(self, v): self.details[-2] = v
+
+    @property
+    def theta_par(self): return self.details[-1]
+
+    def show(self):
+        print("total_pd", self.total_pd)
+        print("num_active", self.num_active)
+        print("pd_par", self.pd_par)
+        print("pd_length", self.pd_length)
+        print("pd_offset", self.pd_offset)
+        print("pd_stride", self.pd_stride)
+        print("par_offsets", self.par_offset)
+        print("num_coord", self.num_coord)
+        print("par_coord", self.par_coord)
+        print("pd_coord", self.pd_coord)
+        print("theta par", self.details[-1])
+
+def mono_details(model_info):
+    details = CoordinationDetails(model_info)
+    # The zero defaults for monodisperse systems are mostly fine
+    details.par_offset[:] = np.arange(2, len(details.par_offset)+2)
+    return details
+
+def poly_details(model_info, weights):
+    #print("weights",weights)
+    weights = weights[2:] # Skip scale and background
+
+    # 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])
+    num_active = np.sum(pd_length>1)
+    if num_active > model_info['parameters'].max_pd:
+        raise ValueError("Too many polydisperse parameters")
+
+    pd_offset = np.cumsum(np.hstack((0, pd_length)))
+    idx = np.argsort(pd_length)[::-1][:num_active]
+    par_length = np.array([max(len(w),1) for w in weights])
+    pd_stride = np.cumprod(np.hstack((1, par_length[idx])))
+    par_offsets = np.cumsum(np.hstack((2, par_length)))
+
+    details = CoordinationDetails(model_info)
+    details.pd_par[:num_active] = idx
+    details.pd_length[:num_active] = pd_length[idx]
+    details.pd_offset[:num_active] = pd_offset[idx]
+    details.pd_stride[:num_active] = pd_stride[:-1]
+    details.par_offset[:] = par_offsets[:-1]
+    details.total_pd = pd_stride[-1]
+    details.num_active = num_active
+    # Without constraints coordinated parameters are just the pd parameters
+    details.par_coord[:num_active] = idx
+    details.pd_coord[:num_active] = 2**np.arange(num_active)
+    details.num_coord = num_active
+    #details.show()
+    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_vector_Iq(model_info):
+    """
+    Define Iq as a vector function if it exists.
+    """
+    Iq = model_info['Iq']
+    if callable(Iq) and not getattr(Iq, 'vectorized', False):
+        def vector_Iq(q, *args):
+            """
+            Vectorized 1D kernel.
+            """
+            return np.array([Iq(qi, *args) for qi in q])
+        model_info['Iq'] = vector_Iq
+
+def create_vector_Iqxy(model_info):
+    """
+    Define Iqxy as a vector function if it exists, or default it from Iq().
+    """
+    Iq, Iqxy = model_info['Iq'], model_info['Iqxy']
+    if callable(Iqxy) and not getattr(Iqxy, 'vectorized', False):
+        def vector_Iqxy(qx, qy, *args):
+            """
+            Vectorized 2D kernel.
+            """
+            return np.array([Iqxy(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)])
+        model_info['Iqxy'] = vector_Iqxy
+    elif callable(Iq):
+        # Iq is vectorized because create_vector_Iq was already called.
+        def default_Iqxy(qx, qy, *args):
+            """
+            Default 2D kernel.
+            """
+            return Iq(np.sqrt(qx**2 + qy**2), *args)
+        model_info['Iqxy'] = default_Iqxy
+
+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.  This also vectorizes
+    any functions that are not already marked as vectorized.
+    """
+    create_vector_Iq(model_info)
+    create_vector_Iqxy(model_info)  # call create_vector_Iq() first
 
 def load_kernel_module(model_name):
@@ -682,9 +710,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
@@ -700,8 +724,8 @@
       *model_info* blocks for the composition objects.  This allows us to
       build complete product and mixture models from just the info.
-
     """
     # 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]
@@ -713,6 +737,6 @@
         filename=abspath(kernel_module.__file__),
         name=name,
-        title=kernel_module.title,
-        description=kernel_module.description,
+        title=getattr(kernel_module, 'title', name+" model"),
+        description=getattr(kernel_module, 'description', 'no description'),
         parameters=parameters,
         composition=None,
@@ -721,4 +745,5 @@
         single=getattr(kernel_module, 'single', True),
         structure_factor=getattr(kernel_module, 'structure_factor', False),
+        profile_axes=getattr(kernel_module, 'profile_axes', ['x','y']),
         variant_info=getattr(kernel_module, 'invariant_info', None),
         demo=getattr(kernel_module, 'demo', None),
@@ -726,8 +751,13 @@
         tests=getattr(kernel_module, 'tests', []),
         )
-    process_parameters(model_info)
+    set_demo(model_info, getattr(kernel_module, 'demo', None))
     # Check for optional functions
-    functions = "ER VR form_volume Iq Iqxy shape sesans".split()
+    functions = "ER VR form_volume Iq Iqxy profile 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
 
@@ -796,4 +826,5 @@
     Program which prints the source produced by the model.
     """
+    import sys
     if len(sys.argv) <= 1:
         print("usage: python -m sasmodels.generate modelname")
Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
+++ sasmodels/kernel_header.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
@@ -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 6e7ff6d026d856e1eb236a05975a2600ba9ba2c9)
+++ sasmodels/kernel_iq.c	(revision 6e7ff6d026d856e1eb236a05975a2600ba9ba2c9)
@@ -0,0 +1,222 @@
+
+/*
+    ##########################################################
+    #                                                        #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #   !!                                              !!   #
+    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
+    #   !!                                              !!   #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #                                                        #
+    ##########################################################
+*/
+
+#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
+#define _PAR_BLOCK_
+
+typedef struct {
+#if MAX_PD > 0
+    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
+#endif // MAX_PD > 0
+    int32_t par_offset[NPARS];  // offset of par value blocks in the value & weight vector
+    int32_t par_coord[NPARS];   // ids of the coordination parameters
+    int32_t pd_coord[NPARS];    // polydispersity coordination bitvector
+    int32_t num_active;         // number of non-trivial pd loops
+    int32_t total_pd;           // total number of voxels in hypercube
+    int32_t num_coord;          // number of coordinated parameters
+    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 *details,
+    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
+
+  // Fill in the initial variables
+  #ifdef USE_OPENMP
+  #pragma omp parallel for
+  #endif
+  for (int k=0; k < NPARS; k++) {
+    pvec[k] = values[details->par_offset[k]];
+  }
+
+  // If it is the first round initialize the result to zero, otherwise
+  // assume that the previous result has been passed back.
+  // Note: doing this even in the monodisperse case in order to handle the
+  // rare case where the model parameters are invalid and zero is returned.
+  // So slightly increased cost for slightly smaller code size.
+  if (pd_start == 0) {
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    for (int i=0; i < nq+1; i++) {
+      result[i] = 0.0;
+    }
+  }
+
+  // Monodisperse computation
+  if (details->num_active == 0) {
+    #ifdef INVALID
+    if (INVALID(local_values)) { return; }
+    #endif
+
+    const double norm = CALL_VOLUME(local_values);
+    const double scale = values[0];
+    const double background = values[1];
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    result[nq] = norm; // Total volume normalization
+    for (int i=0; i < nq; i++) {
+      double scattering = CALL_IQ(q, i, local_values);
+      result[i] = (norm>0. ? scale*scattering/norm + background : background);
+    }
+    return;
+  }
+
+#if MAX_PD > 0
+  //printf("Entering polydispersity from %d to %d\n", pd_start, pd_stop);
+  // Since we are no longer looping over the entire polydispersity hypercube
+  // for each q, we need to track the normalization values between calls.
+  double norm = 0.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 spherical_correction = 1.0;  // cosine correction for latitude variation
+
+  // Location in the polydispersity hypercube, one index per dimension.
+  local int pd_index[MAX_PD];
+
+  // Location of the coordinated parameters in their own sub-cubes.
+  local int offset[NPARS];
+
+  // Trigger the reset behaviour that happens at the end the fast loop
+  // by setting the initial index >= weight vector length.
+  const int fast_length = details->pd_length[0];
+  pd_index[0] = fast_length;
+
+  // 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] == fast_length) {
+      //printf("should be here with %d active\n", details->num_active);
+
+      // Compute position in polydispersity hypercube
+      for (int k=0; k < details->num_active; k++) {
+        pd_index[k] = (loop_index/details->pd_stride[k])%details->pd_length[k];
+        //printf("pd_index[%d] = %d\n",k,pd_index[k]);
+      }
+
+      // Compute partial weights
+      partial_weight = 1.0;
+      //printf("partial weight %d: ", loop_index);
+      for (int k=1; k < details->num_active; k++) {
+        double wi = weights[details->pd_offset[k] + pd_index[k]];
+        //printf("pd[%d]=par[%d]=%g ", k, details->pd_par[k], wi);
+        partial_weight *= wi;
+      }
+      //printf("\n");
+
+      // Update parameter offsets in weight vector
+      //printf("slow %d: ", loop_index);
+      for (int k=0; k < details->num_coord; k++) {
+        int par = details->par_coord[k];
+        int coord = details->pd_coord[k];
+        int this_offset = details->par_offset[par];
+        int block_size = 1;
+        for (int bit=0; coord != 0; bit++) {
+          if (coord&1) {
+              this_offset += block_size * pd_index[bit];
+              block_size *= details->pd_length[bit];
+          }
+          coord >>= 1;
+        }
+        offset[par] = this_offset;
+        pvec[par] = values[this_offset];
+        //printf("par[%d]=v[%d]=%g \n", k, offset[k], pvec[k]);
+        // if theta is not coordinated with fast index, precompute spherical correction
+        if (par == details->theta_par && !(details->par_coord[k]&1)) {
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1e-6);
+        }
+      }
+      //printf("\n");
+    }
+
+    // Increment fast index
+    const double wi = weights[details->pd_offset[0] + pd_index[0]++];
+    double weight = partial_weight*wi;
+    //printf("fast %d: ", loop_index);
+    for (int k=0; k < details->num_coord; k++) {
+      if (details->pd_coord[k]&1) {
+        const int par = details->par_coord[k];
+        pvec[par] = values[offset[par]++];
+        //printf("p[%d]=v[%d]=%g ", par, offset[par]-1, pvec[par]);
+        // if theta is coordinated with fast index, compute spherical correction each time
+        if (par == details->theta_par) {
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1e-6);
+        }
+      }
+    }
+    //printf("\n");
+
+    #ifdef INVALID
+    if (INVALID(local_values)) continue;
+    #endif
+
+    // Accumulate I(q)
+    // Note: weight==0 must always be excluded
+    if (weight > cutoff) {
+      // spherical correction has some nasty effects when theta is +90 or -90
+      // where it becomes zero.  If the entirety of the correction
+      weight *= spherical_correction;
+      norm += weight * 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] += weight*scattering;
+      }
+    }
+  }
+
+  // Accumulate norm.
+  result[nq] += norm;
+
+  // End of the PD loop we can normalize
+  if (pd_stop >= details->total_pd) {
+    const double scale = values[0];
+    const double background = values[1];
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    for (int i=0; i < nq; i++) {
+      result[i] = (norm>0. ? scale*result[i]/norm + background : background);
+    }
+  }
+#endif // MAX_PD > 0
+}
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 4d767110d389aa65268b4b929ac8df1413460380)
+++ sasmodels/kernelcl.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -48,4 +48,5 @@
 harmless, albeit annoying.
 """
+from __future__ import print_function
 import os
 import warnings
@@ -54,4 +55,5 @@
 
 try:
+    raise NotImplementedError("OpenCL not yet implemented for new kernel template")
     import pyopencl as cl
     # Ask OpenCL for the default context so that we know that one exists
@@ -73,4 +75,19 @@
 # of polydisperse parameters.
 MAX_LOOPS = 2048
+
+
+# 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
+"""
 
 
@@ -142,10 +159,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_list)
     program = cl.Program(context, source).build(options=options)
     return program
@@ -220,7 +244,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]
@@ -317,6 +342,6 @@
         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 +390,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 +413,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 +441,14 @@
     """
     def __init__(self, kernel, model_info, q_vectors, dtype):
+        max_pd = model_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 +456,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 4d767110d389aa65268b4b929ac8df1413460380)
+++ sasmodels/kerneldll.py	(revision 4d767110d389aa65268b4b929ac8df1413460380)
@@ -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"
@@ -90,18 +89,18 @@
 
 
-def dll_path(model_info, dtype="double"):
-    """
-    Path to the compiled model defined by *model_info*.
-    """
-    from os.path import join as joinpath, split as splitpath, splitext
-    basename = splitext(splitpath(model_info['filename'])[1])[0]
-    if np.dtype(dtype) == generate.F32:
-        basename += "32"
-    elif np.dtype(dtype) == generate.F64:
-        basename += "64"
-    else:
-        basename += "128"
-    return joinpath(DLL_PATH, basename+'.so')
-
+def dll_name(model_info, dtype):
+    """
+    Name of the dll containing the model.  This is the base file name without
+    any path or extension, with a form such as 'sas_sphere32'.
+    """
+    bits = 8*dtype.itemsize
+    return "sas_%s%d"%(model_info['id'], bits)
+
+def dll_path(model_info, dtype):
+    """
+    Complete path to the dll for the model.  Note that the dll may not
+    exist yet if it hasn't been compiled.
+    """
+    return os.path.join(DLL_PATH, dll_name(model_info, dtype)+".so")
 
 def make_dll(source, model_info, dtype="double"):
@@ -133,18 +132,10 @@
         dtype = generate.F64  # Force 64-bit dll
 
-    if dtype == generate.F32: # 32-bit dll
-        tempfile_prefix = 'sas_' + model_info['name'] + '32_'
-    elif dtype == generate.F64:
-        tempfile_prefix = 'sas_' + model_info['name'] + '64_'
-    else:
-        tempfile_prefix = 'sas_' + model_info['name'] + '128_'
- 
     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
-        fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix)
+        basename = dll_name(model_info, dtype) + "_"
+        fid, filename = tempfile.mkstemp(suffix=".c", prefix=basename)
         os.fdopen(fid, "w").write(source)
         command = COMPILE%{"source":filename, "output":dll}
@@ -171,8 +162,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 +184,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 +194,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):
@@ -272,36 +252,35 @@
     """
     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 isinstance(details, generate.CoordinationDetails)
+        assert weights.dtype == real and values.dtype == real
+
+        start, stop = 0, details.total_pd
+        #print("in kerneldll")
+        #print("weights", weights)
+        #print("values", values)
+        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 4d767110d389aa65268b4b929ac8df1413460380)
+++ sasmodels/kernelpy.py	(revision 6e7ff6d026d856e1eb236a05975a2600ba9ba2c9)
@@ -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):
@@ -82,61 +87,63 @@
     """
     def __init__(self, kernel, model_info, q_input):
+        self.dtype = np.dtype('d')
         self.info = model_info
         self.q_input = q_input
         self.res = np.empty(q_input.nq, q_input.dtype)
-        dim = '2d' if q_input.is_2d else '1d'
-        # Loop over q unless user promises that the kernel is vectorized by
-        # taggining it with vectorized=True
-        if not getattr(kernel, 'vectorized', False):
-            if dim == '2d':
-                def vector_kernel(qx, qy, *args):
-                    """
-                    Vectorized 2D kernel.
-                    """
-                    return np.array([kernel(qxi, qyi, *args)
-                                     for qxi, qyi in zip(qx, qy)])
+        self.kernel = kernel
+        self.dim = '2d' if q_input.is_2d else '1d'
+
+        partable = model_info['parameters']
+        kernel_parameters = (partable.iqxy_parameters if q_input.is_2d
+                             else partable.iq_parameters)
+        volume_parameters = partable.form_volume_parameters
+
+        # Create an array to hold the parameter values.  There will be a
+        # single array whose values are updated as the calculator goes
+        # through the loop.  Arguments to the kernel and volume functions
+        # will use views into this vector, relying on the fact that a
+        # an array of no dimensions acts like a scalar.
+        parameter_vector = np.empty(len(partable.call_parameters), 'd')
+
+        # Create views into the array to hold the arguments
+        offset = 2
+        kernel_args, volume_args = [], []
+        for p in partable.kernel_parameters:
+            if p.length == 1:
+                # Scalar values are length 1 vectors with no dimensions.
+                v = parameter_vector[offset:offset+1].reshape(())
             else:
-                def vector_kernel(q, *args):
-                    """
-                    Vectorized 1D kernel.
-                    """
-                    return np.array([kernel(qi, *args)
-                                     for qi in q])
-            self.kernel = vector_kernel
+                # Vector values are simple views.
+                v = parameter_vector[offset:offset+p.length]
+            offset += p.length
+            if p in kernel_parameters:
+                kernel_args.append(v)
+            if p in volume_parameters:
+                volume_args.append(p)
+
+        # Hold on to the parameter vector so we can use it to call kernel later.
+        # This may also be required to preserve the views into the vector.
+        self._parameter_vector = parameter_vector
+
+        # Generate a closure which calls the kernel with the views into the
+        # parameter array.
+        if q_input.is_2d:
+            form = model_info['Iqxy']
+            qx, qy = q_input.q[:,0], q_input.q[:,1]
+            self._form = lambda: form(qx, qy, *kernel_args)
         else:
-            self.kernel = kernel
-        fixed_pars = model_info['partype']['fixed-' + dim]
-        pd_pars = model_info['partype']['pd-' + dim]
-        vol_pars = model_info['partype']['volume']
-
-        # First two fixed pars are scale and background
-        pars = [p.name for p in model_info['parameters'][2:]]
-        offset = len(self.q_input.q_vectors)
-        self.args = self.q_input.q_vectors + [None] * len(pars)
-        self.fixed_index = np.array([pars.index(p) + offset
-                                     for p in fixed_pars[2:]])
-        self.pd_index = np.array([pars.index(p) + offset
-                                  for p in pd_pars])
-        self.vol_index = np.array([pars.index(p) + offset
-                                   for p in vol_pars])
-        try: self.theta_index = pars.index('theta') + offset
-        except ValueError: self.theta_index = -1
-
-        # Caller needs fixed_pars and pd_pars
-        self.fixed_pars = fixed_pars
-        self.pd_pars = pd_pars
-
-    def __call__(self, fixed, pd, cutoff=1e-5):
-        #print("fixed",fixed)
-        #print("pd", pd)
-        args = self.args[:]  # grab a copy of the args
-        form, form_volume = self.kernel, self.info['form_volume']
-        # First two fixed
-        scale, background = fixed[:2]
-        for index, value in zip(self.fixed_index, fixed[2:]):
-            args[index] = float(value)
-        res = _loops(form, form_volume, cutoff, scale, background, args,
-                     pd, self.pd_index, self.vol_index, self.theta_index)
-
+            form = model_info['Iq']
+            q = q_input.q
+            self._form = lambda: form(q, *kernel_args)
+
+        # Generate a closure which calls the form_volume if it exists.
+        form_volume = model_info['form_volume']
+        self._volume = ((lambda: form_volume(*volume_args)) if form_volume
+                        else (lambda: 1.0))
+
+    def __call__(self, details, weights, values, cutoff):
+        # type: (.generate.CoordinationDetails, np.ndarray, np.ndarray, float) -> np.ndarray
+        res = _loops(self._parameter_vector, self._form, self._volume,
+                     self.q_input.nq, details, weights, values, cutoff)
         return res
 
@@ -147,34 +154,13 @@
         self.q_input = None
 
-def _loops(form, form_volume, cutoff, scale, background,
-           args, pd, pd_index, vol_index, theta_index):
-    """
-    *form* is the name of the form function, which should be vectorized over
-    q, but otherwise have an interface like the opencl kernels, with the
-    q parameters first followed by the individual parameters in the order
-    given in model.parameters (see :mod:`sasmodels.generate`).
-
-    *form_volume* calculates the volume of the shape.  *vol_index* gives
-    the list of volume parameters
-
-    *cutoff* ignores the corners of the dispersion hypercube
-
-    *scale*, *background* multiplies the resulting form and adds an offset
-
-    *args* is the prepopulated set of arguments to the form function, starting
-    with the q vectors, and including slots for all the parameters.  The
-    values for the parameters will be substituted with values from the
-    dispersion functions.
-
-    *pd* is the list of dispersion parameters
-
-    *pd_index* are the indices of the dispersion parameters in *args*
-
-    *vol_index* are the indices of the volume parameters in *args*
-
-    *theta_index* is the index of the theta parameter for the sasview
-    spherical correction, or -1 if there is no angular dispersion
-    """
-
+def _loops(parameters,  # type: np.ndarray
+           form,        # type: Callable[[], np.ndarray]
+           form_volume, # type: Callable[[], float]
+           nq,          # type: int
+           details,     # type: .generate.CoordinationDetails
+           weights,     # type: np.ndarray
+           values,      # type: np.ndarray
+           cutoff,      # type: float
+           ):           # type: (...) -> None
     ################################################################
     #                                                              #
@@ -186,76 +172,70 @@
     #                                                              #
     ################################################################
-
-    weight = np.empty(len(pd), 'd')
-    if weight.size > 0:
-        # weight vector, to be populated by polydispersity loops
-
-        # identify which pd parameters are volume parameters
-        vol_weight_index = np.array([(index in vol_index) for index in pd_index])
-
-        # Sort parameters in decreasing order of pd length
-        Npd = np.array([len(pdi[0]) for pdi in pd], 'i')
-        order = np.argsort(Npd)[::-1]
-        stride = np.cumprod(Npd[order])
-        pd = [pd[index] for index in order]
-        pd_index = pd_index[order]
-        vol_weight_index = vol_weight_index[order]
-
-        fast_value = pd[0][0]
-        fast_weight = pd[0][1]
+    parameters[2:] = values[details.par_offset]
+    scale, background = values[0], values[1]
+    if details.num_active == 0:
+        norm = float(form_volume())
+        if norm > 0.0:
+            return (scale/norm)*form() + background
+        else:
+            return np.ones(nq, 'd')*background
+
+    partial_weight = np.NaN
+    spherical_correction = 1.0
+    pd_stride = details.pd_stride[:details.num_active]
+    pd_length = details.pd_length[:details.num_active]
+    pd_offset = details.pd_offset[:details.num_active]
+    pd_index = np.empty_like(pd_offset)
+    offset = np.empty_like(details.par_offset)
+    theta = details.theta_par
+    fast_length = pd_length[0]
+    pd_index[0] = fast_length
+    total = np.zeros(nq, 'd')
+    norm = 0.0
+    for loop_index in range(details.total_pd):
+        # update polydispersity parameter values
+        if pd_index[0] == fast_length:
+            pd_index[:] = (loop_index/pd_stride)%pd_length
+            partial_weight = np.prod(weights[pd_offset+pd_index][1:])
+            for k in range(details.num_coord):
+                par = details.par_coord[k]
+                coord = details.pd_coord[k]
+                this_offset = details.par_offset[par]
+                block_size = 1
+                for bit in xrange(32):
+                    if coord&1:
+                        this_offset += block_size * pd_index[bit]
+                        block_size *= pd_length[bit]
+                    coord >>= 1
+                    if coord == 0: break
+                offset[par] = this_offset
+                parameters[par] = values[this_offset]
+                if par == theta and not (details.par_coord[k]&1):
+                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6)
+        for k in range(details.num_coord):
+            if details.pd_coord[k]&1:
+                par = details.par_coord[k]
+                parameters[par] = values[offset[par]]
+                offset[par] += 1
+                if par == theta:
+                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6)
+
+        weight = partial_weight * weights[pd_offset[0] + pd_index[0]]
+        pd_index[0] += 1
+        if weight > cutoff:
+            # Call the scattering function
+            # Assume that NaNs are only generated if the parameters are bad;
+            # exclude all q for that NaN.  Even better would be to have an
+            # INVALID expression like the C models, but that is too expensive.
+            I = form()
+            if np.isnan(I).any(): continue
+
+            # update value and norm
+            weight *= spherical_correction
+            total += weight * I
+            norm += weight * form_volume()
+
+    if norm > 0.0:
+        return (scale/norm)*total + background
     else:
-        stride = np.array([1])
-        vol_weight_index = slice(None, None)
-        # keep lint happy
-        fast_value = [None]
-        fast_weight = [None]
-
-    ret = np.zeros_like(args[0])
-    norm = np.zeros_like(ret)
-    vol = np.zeros_like(ret)
-    vol_norm = np.zeros_like(ret)
-    for k in range(stride[-1]):
-        # update polydispersity parameter values
-        fast_index = k % stride[0]
-        if fast_index == 0:  # bottom loop complete ... check all other loops
-            if weight.size > 0:
-                for i, index, in enumerate(k % stride):
-                    args[pd_index[i]] = pd[i][0][index]
-                    weight[i] = pd[i][1][index]
-        else:
-            args[pd_index[0]] = fast_value[fast_index]
-            weight[0] = fast_weight[fast_index]
-
-        # Computes the weight, and if it is not sufficient then ignore this
-        # parameter set.
-        # Note: could precompute w1*...*wn so we only need to multiply by w0
-        w = np.prod(weight)
-        if w > cutoff:
-            # Note: can precompute spherical correction if theta_index is not
-            # the fast index. Correction factor for spherical integration
-            #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0
-            spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2
-                                    if theta_index >= 0 else 1.0)
-            #spherical_correction = 1.0
-
-            # Call the scattering function and adds it to the total.
-            I = form(*args)
-            if np.isnan(I).any(): continue
-            ret += w * I * spherical_correction
-            norm += w
-
-            # Volume normalization.
-            # If there are "volume" polydispersity parameters, then these
-            # will be used to call the form_volume function from the user
-            # supplied kernel, and accumulate a normalized weight.
-            # Note: can precompute volume norm if fast index is not a volume
-            if form_volume:
-                vol_args = [args[index] for index in vol_index]
-                vol_weight = np.prod(weight[vol_weight_index])
-                vol += vol_weight * form_volume(*vol_args)
-                vol_norm += vol_weight
-
-    positive = (vol * vol_norm != 0.0)
-    ret[positive] *= vol_norm[positive] / vol[positive]
-    result = scale * ret / norm + background
-    return result
+        return np.ones(nq, 'd')*background
Index: sasmodels/list_pars.py
===================================================================
--- sasmodels/list_pars.py	(revision a84a0cac83cecece77ff644df4f2552ef1452b62)
+++ sasmodels/list_pars.py	(revision 21b116f73d822594145622e342e0308fc3ecd14a)
@@ -25,5 +25,5 @@
     for name in sorted(MODELS):
         model_info = load_model_info(name)
-        for p in model_info['parameters']:
+        for p in model_info['parameters'].kernel_parameters:
             partable.setdefault(p.name, [])
             partable[p.name].append(name)
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/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision 6e7ff6d026d856e1eb236a05975a2600ba9ba2c9)
+++ sasmodels/modelinfo.py	(revision 6e7ff6d026d856e1eb236a05975a2600ba9ba2c9)
@@ -0,0 +1,518 @@
+
+import numpy as np
+
+# TODO: turn ModelInfo into a proper class
+ModelInfo = dict
+
+MAX_PD = 4
+
+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)
+    return partable
+
+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 = [0, len(choices)-1]
+    else:
+        choices = []
+    # 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 Exception:
+        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 Exception:
+            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
+
+
+def set_demo(model_info, demo):
+    """
+    Assign demo parameters to model_info['demo']
+
+    If demo is not defined, then use defaults.
+
+    If demo is defined, override defaults with value from demo.
+
+    If demo references vector fields, such as thickness[n], then support
+    different ways of assigning the demo values, including assigning a
+    specific value (e.g., thickness3=50.0), assigning a new value to all
+    (e.g., thickness=50.0) or assigning values using list notation.
+    """
+    partable = model_info['parameters']
+    if demo is None:
+        result = partable.defaults
+    else:
+        pars = dict((p.id, p) for p in partable.kernel_parameters)
+        result = partable.defaults.copy()
+        vectors = dict((name,value) for name,value in demo.items()
+                       if name in pars and pars[name].length > 1)
+        if vectors:
+            scalars = dict((name, value) for name, value in demo.items()
+                           if name not in pars or pars[name].length == 1)
+            for name, value in vectors.items():
+                if np.isscalar(value):
+                    # support for the form
+                    #    demo(thickness=0, thickness2=50)
+                    for k in range(1, pars[name].length+1):
+                        key = name+str(k)
+                        if key not in scalars:
+                            scalars[key] = vectors
+                else:
+                    # supoprt for the form
+                    #    demo(thickness=[20,10,3])
+                    for (k,v) in enumerate(value):
+                        scalars[name+str(k)] = v
+            result.update(scalars)
+        else:
+            result.update(demo)
+
+    model_info['demo'] = result
+
+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. The volume parameters are used for calls
+    to form_volume within the kernel (required for volume normalization)
+    and for calls to ER and VR for effective radius and volume ratio
+    respectively.
+
+    *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
+
+    * *is_control* is True if the parameter is a control parameter for a vector
+
+    * *polydisperse* is true if the parameter accepts a polydispersity
+
+    * *relative_pd* is true if that polydispersity is a portion of the
+    value (so a 10% length dipsersity would use a polydispersity value of 0.1)
+    rather than absolute dispersisity (such as an angle plus or minus
+    15 degrees).
+
+    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.split('[')[0].strip()
+        self.name = name
+        self.units = units
+        self.default = default
+        self.limits = limits
+        self.type = type
+        self.description = description
+
+        # Length and length_control will be filled in once the complete
+        # parameter table is available.
+        self.length = 1
+        self.length_control = None
+        self.is_control = False
+
+        # TODO: need better control over whether a parameter is polydisperse
+        self.polydisperse = False
+        self.relative_pd = False
+
+        # choices are also set externally.
+        self.choices = []
+
+    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):
+    """
+    ParameterTable manages the list of available parameters.
+
+    There are a couple of complications which mean that the list of parameters
+    for the kernel differs from the list of parameters that the user sees.
+
+    (1) Common parameters.  Scale and background are implicit to every model,
+    but are not passed to the kernel.
+
+    (2) Vector parameters.  Vector parameters are passed to the kernel as a
+    pointer to an array, e.g., thick[], but they are seen by the user as n
+    separate parameters thick1, thick2, ...
+
+    Therefore, the parameter table is organized by how it is expected to be
+    used. The following information is needed to set up the kernel functions:
+
+    * *kernel_parameters* is the list of parameters in the kernel parameter
+    table, with vector parameter p declared as p[].
+
+    * *iq_parameters* is the list of parameters to the Iq(q, ...) function,
+    with vector parameter p sent as p[].
+
+    * *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...)
+    function, with vector parameter p sent as p[].
+
+    * *form_volume_parameters* is the list of parameters to the form_volume(...)
+    function, with vector parameter p sent as p[].
+
+    Problem details, which sets up the polydispersity loops, requires the
+    following:
+
+    * *theta_offset* is the offset of the theta parameter in the kernel parameter
+    table, with vector parameters counted as n individual parameters
+    p1, p2, ..., or offset is -1 if there is no theta parameter.
+
+    * *max_pd* is the maximum number of polydisperse parameters, with vector
+    parameters counted as n individual parameters p1, p2, ...  Note that
+    this number is limited to sasmodels.modelinfo.MAX_PD.
+
+    * *npars* is the total number of parameters to the kernel, with vector
+    parameters counted as n individual parameters p1, p2, ...
+
+    * *call_parameters* is the complete list of parameters to the kernel,
+    including scale and background, with vector parameters recorded as
+    individual parameters p1, p2, ...
+
+    * *active_1d* is the set of names that may be polydisperse for 1d data
+
+    * *active_2d* is the set of names that may be polydisperse for 2d data
+
+    User parameters are the set of parameters visible to the user, including
+    the scale and background parameters that the kernel does not see.  User
+    parameters don't use vector notation, and instead use p1, p2, ...
+
+    * *control_parameters* is the
+
+    """
+    # scale and background are implicit parameters
+    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS]
+
+    def __init__(self, parameters):
+        self.kernel_parameters = parameters
+        self._set_vector_lengths()
+        self._make_call_parameter_list()
+        self._categorize_parameters()
+        self._set_defaults()
+        #self._name_table= dict((p.id, p) for p in parameters)
+
+    def _set_vector_lengths(self):
+        # Sort out the length of the vector parameters such as thickness[n]
+        for p in self.kernel_parameters:
+            if p.length_control:
+                for ref in self.kernel_parameters:
+                    if ref.id == p.length_control:
+                        break
+                else:
+                    raise ValueError("no reference variable %r for %s"
+                                     % (p.length_control, p.name))
+                ref.is_control = True
+                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 = high
+
+    def _set_defaults(self):
+        # Construct default values, including vector defaults
+        defaults = {}
+        for p in self.call_parameters:
+            if p.length == 1:
+                defaults[p.id] = p.default
+            else:
+                for k in range(1, p.length+1):
+                    defaults["%s%d"%(p.id, k)] = p.default
+        self.defaults = defaults
+
+    def _make_call_parameter_list(self):
+        full_list = self.COMMON[:]
+        for p in self.kernel_parameters:
+            if p.length == 1:
+                full_list.append(p)
+            else:
+                for k in range(1, p.length+1):
+                    pk = Parameter(p.id+str(k), p.units, p.default,
+                                   p.limits, p.type, p.description)
+                    pk.polydisperse = p.polydisperse
+                    pk.relative_pd = p.relative_pd
+                    full_list.append(pk)
+        self.call_parameters = full_list
+
+    def _categorize_parameters(self):
+        # Set the kernel parameters.  Assumes background and scale are the
+        # first two parameters in the parameter list, but these are not sent
+        # to the underlying kernel functions.
+        self.iq_parameters = [p for p in self.kernel_parameters
+                              if p.type not in ('orientation', 'magnetic')]
+        self.iqxy_parameters = [p for p in self.kernel_parameters
+                                if p.type != 'magnetic']
+        self.form_volume_parameters = [p for p in self.kernel_parameters
+                                       if p.type == 'volume']
+
+        # Theta offset
+        offset = 0
+        for p in self.kernel_parameters:
+            if p.name == 'theta':
+                self.theta_offset = offset
+                break
+            offset += p.length
+        else:
+            self.theta_offset = -1
+
+        # number of polydisperse parameters
+        num_pd = sum(p.length for p in self.kernel_parameters if p.polydisperse)
+        # 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.
+        self.max_pd = min(num_pd, MAX_PD)
+
+        self.npars = sum(p.length for p in self.kernel_parameters)
+
+        # true if has 2D parameters
+        self.has_2d = any(p.type in ('orientation', 'magnetic')
+                          for p in self.kernel_parameters)
+
+        self.pd_1d = set(p.name for p in self.call_parameters
+                if p.polydisperse and p.type not in ('orientation', 'magnetic'))
+        self.pd_2d = set(p.name for p in self.call_parameters
+                         if p.polydisperse and p.type != 'magnetic')
+
+    def user_parameters(self, pars={}, is2d=True):
+        """
+        Return the list of parameters for the given data type.
+
+        Vector parameters are expanded as in place.  If multiple parameters
+        share the same vector length, then the parameters will be interleaved
+        in the result.  The control parameters come first.  For example,
+        if the parameter table is ordered as::
+
+            sld_core
+            sld_shell[num_shells]
+            sld_solvent
+            thickness[num_shells]
+            num_shells
+
+        and *pars[num_shells]=2* then the returned list will be::
+
+            num_shells
+            scale
+            background
+            sld_core
+            sld_shell1
+            thickness1
+            sld_shell2
+            thickness2
+            sld_solvent
+
+        Note that shell/thickness pairs are grouped together in the result
+        even though they were not grouped in the incoming table.  The control
+        parameter is always returned first since the GUI will want to set it
+        early, and rerender the table when it is changed.
+        """
+        control = [p for p in self.kernel_parameters if p.is_control]
+
+        # Gather entries such as name[n] into groups of the same n
+        dependent = dict((p.id, []) for p in control)
+        for p in self.kernel_parameters:
+            if p.length_control is not None:
+                dependent[p.length_control].append(p)
+
+        # Gather entries such as name[4] into groups of the same length
+        fixed = {}
+        for p in self.kernel_parameters:
+            if p.length > 1 and p.length_control is None:
+                fixed.setdefault(p.length, []).append(p)
+
+        # Using the call_parameters table, we already have expanded forms
+        # for each of the vector parameters; put them in a lookup table
+        expanded_pars = dict((p.name, p) for p in self.call_parameters)
+
+        # Gather the user parameters in order
+        result = control + self.COMMON
+        for p in self.kernel_parameters:
+            if not is2d and p.type in ('orientation', 'magnetic'):
+                pass
+            elif p.is_control:
+                pass # already added
+            elif p.length_control is not None:
+                table = dependent.get(p.length_control, [])
+                if table:
+                    # look up length from incoming parameters
+                    table_length = int(pars.get(p.length_control, p.length))
+                    del dependent[p.length_control] # first entry seen
+                    for k in range(1, table_length+1):
+                        for entry in table:
+                            result.append(expanded_pars[entry.id+str(k)])
+                else:
+                    pass # already processed all entries
+            elif p.length > 1:
+                table = fixed.get(p.length, [])
+                if table:
+                    table_length = p.length
+                    del fixed[p.length]
+                    for k in range(1, table_length+1):
+                        for entry in table:
+                            result.append(expanded_pars[entry.id+str(k)])
+                else:
+                    pass # already processed all entries
+            else:
+                result.append(p)
+
+        return result
+
+
+
+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/flexible_cylinder.c
===================================================================
--- sasmodels/models/flexible_cylinder.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
+++ sasmodels/models/flexible_cylinder.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
@@ -1,59 +1,20 @@
-double form_volume(double length, double kuhn_length, double radius);
-double Iq(double q, double length, double kuhn_length, double radius,
-          double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double length, double kuhn_length,
-            double radius, double sld, double solvent_sld);
-double flexible_cylinder_kernel(double q, double length, double kuhn_length,
-                                double radius, double sld, double solvent_sld);
-
-
-double form_volume(double length, double kuhn_length, double radius)
+static double
+form_volume(length, kuhn_length, radius)
 {
     return 1.0;
 }
 
-double flexible_cylinder_kernel(double q,
-          double length,
-          double kuhn_length,
-          double radius,
-          double sld,
-          double solvent_sld)
+static double
+Iq(double q,
+   double length,
+   double kuhn_length,
+   double radius,
+   double sld,
+   double solvent_sld)
 {
-
-    const double cont = sld-solvent_sld;
-    const double qr = q*radius;
-    //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr);
-    const double crossSect = sas_J1c(qr);
-    double flex = Sk_WR(q,length,kuhn_length);
-    flex *= crossSect*crossSect;
-    flex *= M_PI*radius*radius*length;
-    flex *= cont*cont;
-    flex *= 1.0e-4;
-    return flex;
+    const double contrast = sld - solvent_sld;
+    const double cross_section = sas_J1c(q*radius);
+    const double volume = M_PI*radius*radius*length;
+    const double flex = Sk_WR(q, length, kuhn_length);
+    return 1.0e-4 * volume * square(contrast*cross_section) * flex;
 }
-
-double Iq(double q,
-          double length,
-          double kuhn_length,
-          double radius,
-          double sld,
-          double solvent_sld)
-{
-
-    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
-    return result;
-}
-
-double Iqxy(double qx, double qy,
-            double length,
-            double kuhn_length,
-            double radius,
-            double sld,
-            double solvent_sld)
-{
-    double q;
-    q = sqrt(qx*qx+qy*qy);
-    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
-
-    return result;
-}
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/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/hardsphere.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -149,9 +149,4 @@
    """
 
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/hayter_msa.py
===================================================================
--- sasmodels/models/hayter_msa.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/hayter_msa.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -87,8 +87,4 @@
     return 1.0;
     """
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar.py
===================================================================
--- sasmodels/models/lamellar.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lamellar.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -82,8 +82,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_hg.py
===================================================================
--- sasmodels/models/lamellar_hg.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lamellar_hg.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -101,8 +101,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_hg_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_hg_stack_caille.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lamellar_hg_stack_caille.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -120,8 +120,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_stack_caille.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lamellar_stack_caille.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -104,8 +104,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_stack_paracrystal.py
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lamellar_stack_paracrystal.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -132,8 +132,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lib/sas_JN.c
===================================================================
--- sasmodels/models/lib/sas_JN.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
+++ sasmodels/models/lib/sas_JN.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
@@ -48,6 +48,7 @@
 */
 
-static double
-sas_JN( int n, double x ) {
+double sas_JN( int n, double x );
+
+double sas_JN( int n, double x ) {
 
     const double MACHEP = 1.11022302462515654042E-16;
Index: sasmodels/models/lib/sph_j1c.c
===================================================================
--- sasmodels/models/lib/sph_j1c.c	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ sasmodels/models/lib/sph_j1c.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
@@ -7,4 +7,5 @@
 * using double precision that are the source.
 */
+double sph_j1c(double q);
 
 // The choice of the number of terms in the series and the cutoff value for
@@ -43,5 +44,4 @@
 #endif
 
-double sph_j1c(double q);
 double sph_j1c(double q)
 {
Index: sasmodels/models/lib/sphere_form.c
===================================================================
--- sasmodels/models/lib/sphere_form.c	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
+++ sasmodels/models/lib/sphere_form.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
@@ -1,10 +1,11 @@
-inline double
-sphere_volume(double radius)
+double sphere_volume(double radius);
+double sphere_form(double q, double radius, double sld, double solvent_sld);
+
+double sphere_volume(double radius)
 {
     return M_4PI_3*cube(radius);
 }
 
-inline double
-sphere_form(double q, double radius, double sld, double solvent_sld)
+double sphere_form(double q, double radius, double sld, double solvent_sld)
 {
     const double fq = sphere_volume(radius) * sph_j1c(q*radius);
Index: sasmodels/models/lib/wrc_cyl.c
===================================================================
--- sasmodels/models/lib/wrc_cyl.c	(revision e7678b2f993854334aca4c45e1c7022b2fdbb509)
+++ sasmodels/models/lib/wrc_cyl.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
@@ -2,4 +2,7 @@
     Functions for WRC implementation of flexible cylinders
 */
+double Sk_WR(double q, double L, double b);
+
+
 static double
 AlphaSquare(double x)
@@ -363,5 +366,4 @@
 }
 
-double Sk_WR(double q, double L, double b);
 double Sk_WR(double q, double L, double b)
 {
Index: sasmodels/models/onion.c
===================================================================
--- sasmodels/models/onion.c	(revision fdb148728a4d4bdcb4cd711d87631f3ee41993fd)
+++ sasmodels/models/onion.c	(revision ce896fde473387374a3fcfc3918ce7a13676b16c)
@@ -4,5 +4,5 @@
     double thickness, double A)
 {
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   const double qr = q * r;
   const double alpha = A * r/thickness;
@@ -19,5 +19,5 @@
     double sinqr, cosqr;
     SINCOS(qr, sinqr, cosqr);
-    fun = -3.0(
+    fun = -3.0*(
             ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq
                 - (alpha*sinqr/qr - cosqr)
@@ -32,5 +32,5 @@
 f_linear(double q, double r, double sld, double slope)
 {
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   const double qr = q * r;
   const double bes = sph_j1c(qr);
@@ -52,5 +52,5 @@
 {
   const double bes = sph_j1c(q * r);
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   return sld * vol * bes;
 }
@@ -64,15 +64,15 @@
     r += thickness[i];
   }
-  return 4.0/3.0 * M_PI * r * r * r;
+  return M_4PI_3*cube(r);
 }
 
 static double
-Iq(double q, double core_sld, double core_radius, double solvent_sld,
-    double n, double in_sld[], double out_sld[], double thickness[],
+Iq(double q, double sld_core, double core_radius, double sld_solvent,
+    double n, double sld_in[], double sld_out[], double thickness[],
     double A[])
 {
   int i;
-  r = core_radius;
-  f = f_constant(q, r, core_sld);
+  double r = core_radius;
+  double f = f_constant(q, r, sld_core);
   for (i=0; i<n; i++){
     const double r0 = r;
@@ -92,6 +92,6 @@
     }
   }
-  f -= f_constant(q, r, solvent_sld);
-  f2 = f * f * 1.0e-4;
+  f -= f_constant(q, r, sld_solvent);
+  const double f2 = f * f * 1.0e-4;
 
   return f2;
Index: sasmodels/models/onion.py
===================================================================
--- sasmodels/models/onion.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/onion.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
@@ -293,15 +293,15 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
                "Core scattering length density"],
               ["core_radius", "Ang", 200., [0, inf], "volume",
                "Radius of the core"],
-              ["solvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "",
+              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
                "Solvent scattering length density"],
               ["n", "", 1, [0, 10], "volume",
                "number of shells"],
-              ["in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
+              ["sld_in[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
                "scattering length density at the inner radius of shell k"],
-              ["out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "",
+              ["sld_out[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "",
                "scattering length density at the outer radius of shell k"],
               ["thickness[n]", "Ang", 40., [0, inf], "volume",
@@ -311,14 +311,11 @@
               ]
 
-#source = ["lib/sph_j1c.c", "onion.c"]
-
-def Iq(q, *args, **kw):
-    return q
-
-def Iqxy(qx, *args, **kw):
-    return qx
-
-
-def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A):
+source = ["lib/sph_j1c.c", "onion.c"]
+
+#def Iq(q, *args, **kw):
+#    return q
+
+profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
+def profile(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A):
     """
     SLD profile
@@ -374,10 +371,10 @@
 
 demo = {
-    "solvent_sld": 2.2,
-    "core_sld": 1.0,
+    "sld_solvent": 2.2,
+    "sld_core": 1.0,
     "core_radius": 100,
     "n": 4,
-    "in_sld": [0.5, 1.5, 0.9, 2.0],
-    "out_sld": [nan, 0.9, 1.2, 1.6],
+    "sld_in": [0.5, 1.5, 0.9, 2.0],
+    "sld_out": [nan, 0.9, 1.2, 1.6],
     "thickness": [50, 75, 150, 75],
     "A": [0, -1, 1e-4, 1],
Index: sasmodels/models/rpa.c
===================================================================
--- sasmodels/models/rpa.c	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
+++ sasmodels/models/rpa.c	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -1,32 +1,10 @@
 double Iq(double q, double case_num,
-    double Na, double Phia, double va, double a_sld, double ba,
-    double Nb, double Phib, double vb, double b_sld, double bb,
-    double Nc, double Phic, double vc, double c_sld, double bc,
-    double Nd, double Phid, double vd, double d_sld, double bd,
+    double N[], double Phi[], double v[], double L[], double b[],
     double Kab, double Kac, double Kad,
     double Kbc, double Kbd, double Kcd
     );
 
-double Iqxy(double qx, double qy, double case_num,
-    double Na, double Phia, double va, double a_sld, double ba,
-    double Nb, double Phib, double vb, double b_sld, double bb,
-    double Nc, double Phic, double vc, double c_sld, double bc,
-    double Nd, double Phid, double vd, double d_sld, double bd,
-    double Kab, double Kac, double Kad,
-    double Kbc, double Kbd, double Kcd
-    );
-
-double form_volume(void);
-
-double form_volume(void)
-{
-    return 1.0;
-}
-
 double Iq(double q, double case_num,
-    double Na, double Phia, double va, double La, double ba,
-    double Nb, double Phib, double vb, double Lb, double bb,
-    double Nc, double Phic, double vc, double Lc, double bc,
-    double Nd, double Phid, double vd, double Ld, double bd,
+    double N[], double Phi[], double v[], double L[], double b[],
     double Kab, double Kac, double Kad,
     double Kbc, double Kbd, double Kcd
@@ -36,39 +14,39 @@
 #if 0  // Sasview defaults
   if (icase <= 1) {
-    Na=Nb=1000.0;
-    Phia=Phib=0.0000001;
+    N[0]=N[1]=1000.0;
+    Phi[0]=Phi[1]=0.0000001;
     Kab=Kac=Kad=Kbc=Kbd=-0.0004;
-    La=Lb=1e-12;
-    va=vb=100.0;
-    ba=bb=5.0;
+    L[0]=L[1]=1e-12;
+    v[0]=v[1]=100.0;
+    b[0]=b[1]=5.0;
   } else if (icase <= 4) {
-    Phia=0.0000001;
+    Phi[0]=0.0000001;
     Kab=Kac=Kad=-0.0004;
-    La=1e-12;
-    va=100.0;
-    ba=5.0;
+    L[0]=1e-12;
+    v[0]=100.0;
+    b[0]=5.0;
   }
 #else
   if (icase <= 1) {
-    Na=Nb=0.0;
-    Phia=Phib=0.0;
+    N[0]=N[1]=0.0;
+    Phi[0]=Phi[1]=0.0;
     Kab=Kac=Kad=Kbc=Kbd=0.0;
-    La=Lb=Ld;
-    va=vb=vd;
-    ba=bb=0.0;
+    L[0]=L[1]=L[3];
+    v[0]=v[1]=v[3];
+    b[0]=b[1]=0.0;
   } else if (icase <= 4) {
-    Na = 0.0;
-    Phia=0.0;
+    N[0] = 0.0;
+    Phi[0]=0.0;
     Kab=Kac=Kad=0.0;
-    La=Ld;
-    va=vd;
-    ba=0.0;
+    L[0]=L[3];
+    v[0]=v[3];
+    b[0]=0.0;
   }
 #endif
 
-  const double Xa = q*q*ba*ba*Na/6.0;
-  const double Xb = q*q*bb*bb*Nb/6.0;
-  const double Xc = q*q*bc*bc*Nc/6.0;
-  const double Xd = q*q*bd*bd*Nd/6.0;
+  const double Xa = q*q*b[0]*b[0]*N[0]/6.0;
+  const double Xb = q*q*b[1]*b[1]*N[1]/6.0;
+  const double Xc = q*q*b[2]*b[2]*N[2]/6.0;
+  const double Xd = q*q*b[3]*b[3]*N[3]/6.0;
 
   // limit as Xa goes to 0 is 1
@@ -98,33 +76,33 @@
 #if 0
   const double S0aa = icase<5
-                      ? 1.0 : Na*Phia*va*Paa;
+                      ? 1.0 : N[0]*Phi[0]*v[0]*Paa;
   const double S0bb = icase<2
-                      ? 1.0 : Nb*Phib*vb*Pbb;
-  const double S0cc = Nc*Phic*vc*Pcc;
-  const double S0dd = Nd*Phid*vd*Pdd;
+                      ? 1.0 : N[1]*Phi[1]*v[1]*Pbb;
+  const double S0cc = N[2]*Phi[2]*v[2]*Pcc;
+  const double S0dd = N[3]*Phi[3]*v[3]*Pdd;
   const double S0ab = icase<8
-                      ? 0.0 : sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb;
+                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
   const double S0ac = icase<9
-                      ? 0.0 : sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);
+                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
   const double S0ad = icase<9
-                      ? 0.0 : sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);
+                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
   const double S0bc = (icase!=4 && icase!=7 && icase!= 9)
-                      ? 0.0 : sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc;
+                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
   const double S0bd = (icase!=4 && icase!=7 && icase!= 9)
-                      ? 0.0 : sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);
+                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
   const double S0cd = (icase==0 || icase==2 || icase==5)
-                      ? 0.0 : sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd;
+                      ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
 #else  // sasview equivalent
-//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,Nc,Phic,vc,Pcc);
-  double S0aa = Na*Phia*va*Paa;
-  double S0bb = Nb*Phib*vb*Pbb;
-  double S0cc = Nc*Phic*vc*Pcc;
-  double S0dd = Nd*Phid*vd*Pdd;
-  double S0ab = sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb;
-  double S0ac = sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);
-  double S0ad = sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);
-  double S0bc = sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc;
-  double S0bd = sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);
-  double S0cd = sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd;
+//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc);
+  double S0aa = N[0]*Phi[0]*v[0]*Paa;
+  double S0bb = N[1]*Phi[1]*v[1]*Pbb;
+  double S0cc = N[2]*Phi[2]*v[2]*Pcc;
+  double S0dd = N[3]*Phi[3]*v[3]*Pdd;
+  double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
+  double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
+  double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
+  double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
+  double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
+  double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
 switch(icase){
   case 0:
@@ -311,7 +289,7 @@
   // Note: 1e-13 to convert from fm to cm for scattering length
   const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13;
-  const double Lad = icase<5 ? 0.0 : (La/va - Ld/vd)*sqrt_Nav;
-  const double Lbd = icase<2 ? 0.0 : (Lb/vb - Ld/vd)*sqrt_Nav;
-  const double Lcd = (Lc/vc - Ld/vd)*sqrt_Nav;
+  const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav;
+  const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav;
+  const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav;
 
   const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33
@@ -321,23 +299,2 @@
 
 }
-
-double Iqxy(double qx, double qy,
-    double case_num,
-    double Na, double Phia, double va, double a_sld, double ba,
-    double Nb, double Phib, double vb, double b_sld, double bb,
-    double Nc, double Phic, double vc, double c_sld, double bc,
-    double Nd, double Phid, double vd, double d_sld, double bd,
-    double Kab, double Kac, double Kad,
-    double Kbc, double Kbd, double Kcd
-    )
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return Iq(q,
-        case_num,
-        Na, Phia, va, a_sld, ba,
-        Nb, Phib, vb, b_sld, bb,
-        Nc, Phic, vc, c_sld, bc,
-        Nd, Phid, vd, d_sld, bd,
-        Kab, Kac, Kad,
-        Kbc, Kbd, Kcd);
-}
Index: sasmodels/models/rpa.py
===================================================================
--- sasmodels/models/rpa.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/rpa.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
@@ -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/models/spherical_sld.py
===================================================================
--- sasmodels/models/spherical_sld.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/spherical_sld.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -170,5 +170,5 @@
 # pylint: disable=bad-whitespace, line-too-long
 #            ["name", "units", default, [lower, upper], "type", "description"],
-parameters = [["n_shells",         "",               1,      [0, 9],         "", "number of shells"],
+parameters = [["n",                "",               1,      [0, 9],         "", "number of shells"],
               ["radius_core",      "Ang",            50.0,   [0, inf],       "", "intern layer thickness"],
               ["sld_core",         "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"],
@@ -192,5 +192,5 @@
 
 demo = dict(
-    n_shells=4,
+    n=4,
     scale=1.0,
     solvent_sld=1.0,
Index: sasmodels/models/squarewell.py
===================================================================
--- sasmodels/models/squarewell.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/squarewell.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -123,8 +123,4 @@
 """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/stickyhardsphere.py
===================================================================
--- sasmodels/models/stickyhardsphere.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/stickyhardsphere.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -171,8 +171,4 @@
 """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision f2473145a6baa12aa3ab15421a80d2f5d004935c)
+++ sasmodels/product.py	(revision f2473145a6baa12aa3ab15421a80d2f5d004935c)
@@ -14,5 +14,4 @@
 
 from .core import call_ER_VR
-from .generate import process_parameters
 
 SCALE=0
@@ -58,5 +57,4 @@
     # Iq, Iqxy, form_volume, ER, VR and sesans
     model_info['composition'] = ('product', [p_info, s_info])
-    process_parameters(model_info)
     return model_info
 
@@ -93,6 +91,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 4d767110d389aa65268b4b929ac8df1413460380)
+++ sasmodels/resolution.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -502,9 +502,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
-        raise ValueError("bad parameters: [%s] not in [%s]"%
-                         (", ".join(sorted(extra)), ", ".join(sorted(keys))))
+    par_set = set([p.name for p in form.info['parameters'].call_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(pars.keys()))))
 
     if np.isscalar(width):
@@ -556,9 +557,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
-        raise ValueError("bad parameters: [%s] not in [%s]"%
-                         (", ".join(sorted(extra)), ", ".join(sorted(keys))))
+    par_set = set([p.name for p in form.info['parameters'].call_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(pars.keys()))))
 
     _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 4d767110d389aa65268b4b929ac8df1413460380)
+++ sasmodels/sasview_model.py	(revision ee8f73486ed89c512902019a13b040e576091f49)
@@ -26,4 +26,5 @@
 from . import custom
 from . import generate
+from . import weights
 
 def load_standard_models():
@@ -38,5 +39,5 @@
         try:
             models.append(_make_standard_model(name))
-        except:
+        except Exception:
             logging.error(traceback.format_exc())
     return models
@@ -84,10 +85,19 @@
         self._model = None
         model_info = self._model_info
+        parameters = model_info['parameters']
 
         self.name = model_info['name']
         self.description = model_info['description']
         self.category = None
-        self.multiplicity_info = None
-        self.is_multifunc = False
+        #self.is_multifunc = False
+        for p in parameters.kernel_parameters:
+            if p.is_control:
+                profile_axes = model_info['profile_axes']
+                self.multiplicity_info = [
+                    p.limits[1], p.name, p.choices, profile_axes[0]
+                    ]
+                break
+        else:
+            self.multiplicity_info = []
 
         ## interpret the parameters
@@ -96,24 +106,27 @@
         self.params = collections.OrderedDict()
         self.dispersion = dict()
-        partype = model_info['partype']
-
-        for p in model_info['parameters']:
+
+        self.orientation_params = []
+        self.magnetic_params = []
+        self.fixed = []
+        for p in parameters.user_parameters():
             self.params[p.name] = p.default
             self.details[p.name] = [p.units] + p.limits
-
-        for name in partype['pd-2d']:
-            self.dispersion[name] = {
-                'width': 0,
-                'npts': 35,
-                'nsigmas': 3,
-                'type': 'gaussian',
-            }
-
-        self.orientation_params = (
-            partype['orientation']
-            + [n + '.width' for n in partype['orientation']]
-            + partype['magnetic'])
-        self.magnetic_params = partype['magnetic']
-        self.fixed = [n + '.width' for n in partype['pd-2d']]
+            if p.polydisperse:
+                self.dispersion[p.name] = {
+                    'width': 0,
+                    'npts': 35,
+                    'nsigmas': 3,
+                    'type': 'gaussian',
+                }
+            if p.type == 'orientation':
+                self.orientation_params.append(p.name)
+                self.orientation_params.append(p.name+".width")
+                self.fixed.append(p.name+".width")
+            if p.type == 'magnetic':
+                self.orientation_params.append(p.name)
+                self.magnetic_params.append(p.name)
+                self.fixed.append(p.name+".width")
+
         self.non_fittable = []
 
@@ -234,7 +247,8 @@
         """
         # TODO: fix test so that parameter order doesn't matter
-        ret = ['%s.%s' % (d.lower(), p)
-               for d in self._model_info['partype']['pd-2d']
-               for p in ('npts', 'nsigmas', 'width')]
+        ret = ['%s.%s' % (p.name.lower(), ext)
+               for p in self._model_info['parameters'].user_parameters()
+               for ext in ('npts', 'nsigmas', 'width')
+               if p.polydisperse]
         #print(ret)
         return ret
@@ -309,6 +323,5 @@
             # Check whether we have a list of ndarrays [qx,qy]
             qx, qy = qdist
-            partype = self._model_info['partype']
-            if not partype['orientation'] and not partype['magnetic']:
+            if not self._model_info['parameters'].has_2d:
                 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2))
             else:
@@ -335,10 +348,11 @@
             self._model = core.build_model(self._model_info)
         q_vectors = [np.asarray(q) for q in args]
-        fn = self._model.make_kernel(q_vectors)
-        pars = [self.params[v] for v in fn.fixed_pars]
-        pd_pars = [self._get_weights(p) for p in fn.pd_pars]
-        result = fn(pars, pd_pars, self.cutoff)
-        fn.q_input.release()
-        fn.release()
+        kernel = self._model.make_kernel(q_vectors)
+        pairs = [self._get_weights(p)
+                 for p in self._model_info['parameters'].call_parameters]
+        details, weights, values = core.build_details(kernel, pairs)
+        result = kernel(details, weights, values, cutoff=self.cutoff)
+        kernel.q_input.release()
+        kernel.release()
         return result
 
@@ -415,13 +429,12 @@
         Return dispersion weights for parameter
         """
-        from . import weights
-        relative = self._model_info['partype']['pd-rel']
-        limits = self._model_info['limits']
-        dis = self.dispersion[par]
-        value, weight = weights.get_weights(
-            dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
-            self.params[par], limits[par], par in relative)
-        return value, weight / np.sum(weight)
-
+        if par.polydisperse:
+            dis = self.dispersion[par.name]
+            value, weight = weights.get_weights(
+                dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
+                self.params[par.name], par.limits, par.relative_pd)
+            return value, weight / np.sum(weight)
+        else:
+            return [self.params[par.name]], []
 
 def test_model():
