Index: CylinderModel.py
===================================================================
--- CylinderModel.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
+++ CylinderModel.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -0,0 +1,7 @@
+"""
+Drop-in replacement for sasview cylinder model.
+
+No rescaling or renaming of the parameters.
+"""
+from sasmodels.sasview_model import make_class
+CylinderModel = make_class('cylinder_clone','CylinderModel')
Index: compare-new.py
===================================================================
--- compare-new.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ compare-new.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -48,5 +48,5 @@
             gpu = model.theory()
         gpu_time = toc()*1000./Ngpu
-        print "ocl t=%.1f ms"%gpu_time
+        print "ocl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[~np.isnan(gpu)]))
         #print max(gpu), min(gpu)
 
@@ -61,5 +61,5 @@
         print "dll t=%.1f ms"%cpu_time
 
-    elif 1: # Hack to check new vs old for GpuCylinder
+    elif 0: # Hack to check new vs old for GpuCylinder
         from Models.code_cylinder_f import GpuCylinder as oldgpu
         from sasmodel import SasModel
@@ -78,5 +78,5 @@
             cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data])
         cpu_time = toc()*1000./Ncpu
-        print "sasview t=%.1f ms"%cpu_time
+        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index]))
 
     if Ngpu > 0 and Ncpu > 0:
Index: sasmodel.py
===================================================================
--- sasmodel.py	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
+++ sasmodel.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -4,8 +4,4 @@
 import datetime
 import numpy as np
-import pyopencl as cl
-from bumps.names import Parameter
-from sans.dataloader.loader import Loader
-from sans.dataloader.manipulations import Ringcut, Boxcut
 
 
@@ -20,4 +16,5 @@
 
 def load_data(filename):
+    from sans.dataloader.loader import Loader
     loader = Loader()
     data = loader.load(filename)
@@ -57,4 +54,5 @@
 
 def set_beam_stop(data, radius, outer=None):
+    from sans.dataloader.manipulations import Ringcut
     if hasattr(data, 'qx_data'):
         data.mask = Ringcut(0, radius)(data)
@@ -67,4 +65,5 @@
 
 def set_half(data, half):
+    from sans.dataloader.manipulations import Boxcut
     if half == 'right':
         data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
@@ -73,4 +72,5 @@
 
 def set_top(data, max):
+    from sans.dataloader.manipulations import Boxcut
     data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data)
 
@@ -144,4 +144,5 @@
     global GPU_CONTEXT, GPU_QUEUE
     if GPU_CONTEXT is None:
+        import pyopencl as cl
         GPU_CONTEXT = cl.create_some_context()
         GPU_QUEUE = cl.CommandQueue(GPU_CONTEXT)
@@ -151,4 +152,6 @@
 class SasModel(object):
     def __init__(self, data, model, dtype='float32', **kw):
+        from bumps.names import Parameter
+
         self.__dict__['_parameters'] = {}
         #self.name = data.filename
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/core.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -4,4 +4,5 @@
 import sys, os
 import datetime
+import warnings
 
 import numpy as np
@@ -14,12 +15,14 @@
     return gen.make(modelpath)
 
+
+
 def opencl_model(modelname, dtype="single"):
     from sasmodels import gpu
 
-    source, meta, _ = load_model(modelname)
+    source, info, _ = load_model(modelname)
     # for debugging, save source to a .cl file, edit it, and reload as model
     #open(modelname+'.cl','w').write(source)
     #source = open(modelname+'.cl','r').read()
-    return gpu.GpuModel(source, meta, dtype)
+    return gpu.GpuModel(source, info, dtype)
 
 
@@ -31,8 +34,11 @@
     COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm"
 DLL_PATH = "/tmp"
-def dll_path(meta):
+
+
+def dll_path(info):
     from os.path import join as joinpath, split as splitpath, splitext
-    basename = splitext(splitpath(meta['filename'])[1])[0]
+    basename = splitext(splitpath(info['filename'])[1])[0]
     return joinpath(DLL_PATH, basename+'.so')
+
 
 def dll_model(modelname):
@@ -40,13 +46,14 @@
     from sasmodels import dll
 
-    source, meta, _ = load_model(modelname)
-    dllpath = dll_path(meta)
+    source, info, _ = load_model(modelname)
+    dllpath = dll_path(info)
     if not os.path.exists(dllpath) \
-            or (os.path.getmtime(dllpath) < os.path.getmtime(meta['filename'])):
+            or (os.path.getmtime(dllpath) < os.path.getmtime(info['filename'])):
         # Replace with a proper temp file
         srcfile = '/tmp/%s.c'%modelname
         open(srcfile, 'w').write(source)
         os.system(COMPILE%(srcfile, dllpath))
-    return dll.DllModel(dllpath, meta)
+    return dll.DllModel(dllpath, info)
+
 
 TIC = None
@@ -57,6 +64,8 @@
     return TIC
 
+
 def toc():
     return TIC()
+
 
 def load_data(filename):
@@ -68,7 +77,7 @@
     return data
 
+
 def fake_data2D(qx, qy=None):
     from sans.dataloader.data_info import Data2D, Detector
-
 
     if qy is None:
@@ -121,4 +130,5 @@
             data.mask &= (data.x<outer)
 
+
 def set_half(data, half):
     from sans.dataloader.manipulations import Boxcut
@@ -128,7 +138,9 @@
         data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
 
+
 def set_top(data, max):
     from sans.dataloader.manipulations import Boxcut
     data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data)
+
 
 def plot_data(data, iq, vmin=None, vmax=None):
@@ -141,4 +153,5 @@
                interpolation='nearest', aspect=1, origin='upper',
                extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
+
 
 def plot_result2D(data, theory, view='linear'):
@@ -184,9 +197,9 @@
     mresid = masked_array((theory-data.y)/data.dy, mdata.mask)
 
-    plt.subplot(1,2,1)
+    plt.subplot(121)
     plt.errorbar(data.x, mdata, yerr=data.dy)
     plt.plot(data.x, mtheory, '-', hold=True)
     plt.yscale(view)
-    plt.subplot(1, 2, 2)
+    plt.subplot(122)
     plt.plot(data.x, mresid, 'x')
     #plt.axhline(1, color='black', ls='--',lw=1, hold=True)
@@ -219,26 +232,31 @@
 
         # create model
-        self.fn = model(input, cutoff=cutoff)
+        self.fn = model(input)
+        self.cutoff = cutoff
 
         # define bumps parameters
         pars = []
-        for p in model.meta['parameters']:
+        extras = []
+        for p in model.info['parameters']:
             name, default, limits, ptype = p[0], p[2], p[3], p[4]
             value = kw.pop(name, default)
             setattr(self, name, Parameter.default(value, name=name, limits=limits))
             pars.append(name)
-            if ptype != "":
-                for xpart,xdefault,xlimits in [
-                        ('_pd', 0, limits),
-                        ('_pd_n', 35, (0,1000)),
-                        ('_pd_nsigma', 3, (0,10)),
-                        ]:
-                    xname = name+xpart
-                    xvalue = kw.pop(xname, xdefault)
-                    setattr(self, xname, Parameter.default(xvalue, name=xname))
+        for name in model.info['partype']['pd-2d']:
+            for xpart,xdefault,xlimits in [
+                    ('_pd', 0, limits),
+                    ('_pd_n', 35, (0,1000)),
+                    ('_pd_nsigma', 3, (0, 10)),
+                    ('_pd_type', 'gaussian', None),
+                ]:
+                xname = name+xpart
+                xvalue = kw.pop(xname, xdefault)
+                if xlimits is not None:
+                    xvalue = Parameter.default(xvalue, name=xname, limits=xlimits)
                     pars.append(xname)
+                setattr(self, xname, xvalue)
+        self._parameter_names = pars
         if kw:
             raise TypeError("unexpected parameters: %s"%(", ".join(sorted(kw.keys()))))
-        self._parameter_names = pars
         self.update()
 
@@ -254,8 +272,9 @@
     def theory(self):
         if 'theory' not in self._cache:
-            pars = dict((k,getattr(self,k).value) for k in self._parameter_names)
+            pars = [getattr(self,p).value for p in self.fn.fixed_pars]
+            pd_pars = [self._get_weights(p) for p in self.fn.pd_pars]
             #print pars
-            self._theory[self.index] = self.fn.eval(pars)
-            #self._theory[:] = self.fn.eval(pars)
+            self._theory[self.index] = self.fn(pars, pd_pars, self.cutoff)
+            #self._theory[:] = self.fn.eval(pars, pd_pars)
             self._cache['theory'] = self._theory
         return self._cache['theory']
@@ -282,4 +301,17 @@
         pass
 
+    def _get_weights(self, par):
+        from . import weights
+
+        relative = self.fn.info['partype']['pd-rel']
+        limits = self.fn.info['limits']
+        disperser,value,npts,width,nsigma = [getattr(self, par+ext)
+                for ext in ('_pd_type','','_pd_n','_pd','_pd_nsigma')]
+        v,w = weights.get_weights(
+            disperser, int(npts.value), width.value, nsigma.value,
+            value.value, limits[par], par in relative)
+        return v,w/w.max()
+
+
 def demo():
     data = load_data('DEC07086.DAT')
@@ -288,4 +320,5 @@
     import matplotlib.pyplot as plt; plt.show()
 
+
 if __name__ == "__main__":
     demo()
Index: sasmodels/dll.py
===================================================================
--- sasmodels/dll.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/dll.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -19,5 +19,5 @@
     ctypes wrapper for a single model.
 
-    *source* and *meta* are the model source and interface as returned
+    *source* and *info* are the model source and interface as returned
     from :func:`gen.make`.
 
@@ -27,24 +27,34 @@
     is an optional extension which may not be available on all devices.
     """
-    def __init__(self, dllpath, meta):
-        self.meta = meta
-        self.dll = ct.CDLL(dllpath)
-        self.Iq = self.dll[gen.kernel_name(self.meta, False)]
-        self.Iqxy = self.dll[gen.kernel_name(self.meta, True)]
+    def __init__(self, dllpath, info):
+        self.info = info
+        self.dllpath = dllpath
+        self.dll = None
 
+    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'])
 
-        self.PARS = dict((p[0],p[2]) for p in meta['parameters'])
-        self.PD_PARS = [p[0] for p in meta['parameters'] if p[4] != ""]
+        self.dll = ct.CDLL(self.dllpath)
 
-        # Determine the set of fixed and polydisperse parameters
-        Nfixed = len([p[0] for p in meta['parameters'] if p[4] == ""])
-        N1D = len([p for p in meta['parameters'] if p[4]=="volume"])
-        N2D = len([p for p in meta['parameters'] if p[4]!=""])
-        self.Iq.argtypes = IQ_ARGS + [c_double]*Nfixed + [c_int]*N1D
-        self.Iqxy.argtypes = IQXY_ARGS + [c_double]*Nfixed + [c_int]*N2D
+        self.Iq = self.dll[gen.kernel_name(self.info, False)]
+        self.Iq.argtypes = IQ_ARGS + [c_double]*Nfixed1d + [c_int]*Npd1d
 
-    def __call__(self, input, cutoff=1e-5):
+        self.Iqxy = self.dll[gen.kernel_name(self.info, True)]
+        self.Iqxy.argtypes = IQXY_ARGS + [c_double]*Nfixed2d + [c_int]*Npd2d
+
+    def __getstate__(self):
+        return {'info': self.info, 'dllpath': self.dllpath, 'dll': None}
+
+    def __setstate__(self, state):
+        self.__dict__ = state
+
+    def __call__(self, input):
+        if self.dll is None: self._load_dll()
+
         kernel = self.Iqxy if input.is_2D else self.Iq
-        return DllKernel(kernel, self.meta, input, cutoff)
+        return DllKernel(kernel, self.info, input)
 
     def make_input(self, q_vectors):
@@ -90,29 +100,29 @@
 
 class DllKernel(object):
-    def __init__(self, kernel, meta, input, cutoff):
-        self.cutoff = cutoff
+    def __init__(self, kernel, info, input):
         self.input = input
         self.kernel = kernel
-        self.meta = meta
+        self.info = info
+        self.res = np.empty(input.nq, input.dtype)
+        dim = '2d' if input.is_2D else '1d'
+        self.fixed_pars = info['partype']['fixed-'+dim]
+        self.pd_pars = info['partype']['pd-'+dim]
 
-        self.res = np.empty(input.nq, input.dtype)
+        # In dll kernel, but not in opencl kernel
         self.p_res = self.res.ctypes.data
 
-        # Determine the set of fixed and polydisperse parameters
-        self.fixed_pars = [p[0] for p in meta['parameters'] if p[4] == ""]
-        self.pd_pars = [p for p in meta['parameters']
-               if p[4]=="volume" or (p[4]=="orientation" and input.is_2D)]
+    def __call__(self, pars, pd_pars, cutoff):
+        real = np.float32 if self.input.dtype == F32 else np.float64
+        fixed = [real(p) for p in pars]
+        cutoff = real(cutoff)
+        loops = np.hstack(pd_pars)
+        loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten()
+        loops_N = [np.uint32(len(p[0])) for p in pd_pars]
 
-    def eval(self, pars):
-        fixed, loops, loop_n = \
-            gen.kernel_pars(pars, self.meta, self.input.is_2D, dtype=self.input.dtype)
-        real = np.float32 if self.input.dtype == F32 else np.float64
         nq = c_int(self.input.nq)
-        cutoff = real(self.cutoff)
-
         p_loops = loops.ctypes.data
-        pars = self.input.q_pointers + [self.p_res, nq, p_loops, cutoff] + fixed + loop_n
+        args = self.input.q_pointers + [self.p_res, nq, p_loops, cutoff] + fixed + loops_N
         #print pars
-        self.kernel(*pars)
+        self.kernel(*args)
 
         return self.res
@@ -120,5 +130,2 @@
     def release(self):
         pass
-
-    def __del__(self):
-        self.release()
Index: sasmodels/gen.py
===================================================================
--- sasmodels/gen.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/gen.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -9,4 +9,11 @@
 F64 = np.dtype('float64')
 F32 = np.dtype('float32')
+
+# Scale and background, which are parameters common to every form factor
+COMMON_PARAMETERS = [
+    [ "scale", "", 1, [0, np.inf], "", "Source intensity" ],
+    [ "background", "1/cm", 0, [0, np.inf], "", "Source background" ],
+    ]
+
 
 # Conversion from units defined in the parameter table for each model
@@ -29,11 +36,4 @@
 
 PARTABLE_VALUE_WIDTH = 10
-
-# Scale and background, which are parameters common to every form factor
-COMMON_PARAMETERS = [
-    [ "scale", "", 0, [0, np.inf], "", "Source intensity" ],
-    [ "background", "1/cm", 0, [0, np.inf], "", "Source background" ],
-    ]
-
 
 # Header included before every kernel.
@@ -60,4 +60,5 @@
 #  define kernel
 #  define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0)
+#  define powr(a,b) pow(a,b)
 #else
 #  ifdef USE_SINCOS
@@ -93,5 +94,5 @@
 # respectively, so the template builder will need to do extra work to
 # declare, initialize and pass the q parameters.
-IQ_KERNEL = {
+KERNEL_1D = {
     'fn': "Iq",
     'q_par_decl': "global const real *q,",
@@ -100,5 +101,5 @@
     }
 
-IQXY_KERNEL = {
+KERNEL_2D = {
     'fn': "Iqxy",
     'q_par_decl': "global const real *qx,\n    global const real *qy,",
@@ -176,9 +177,10 @@
 # Volume normalization.
 # If there are "volume" polydispersity parameters, then these will be used
-# to call the volume function from the user supplied kernel, and accumulate
+# to call the form_volume function from the user supplied kernel, and accumulate
 # a normalized weight.
 VOLUME_NORM="""const real vol_weight = %(weight)s;
-  vol += vol_weight*volume(%(pars)s);
+  vol += vol_weight*form_volume(%(pars)s);
   norm_vol += vol_weight;"""
+
 
 def indent(s, depth):
@@ -191,9 +193,13 @@
 
 
-def make_kernel(meta, form):
+def kernel_name(info, is_2D):
+    return info['name'] + "_" + ("Iqxy" if is_2D else "Iq")
+
+
+def make_kernel(info, is_2D):
     """
     Build a kernel call from metadata supplied by the user.
 
-    *meta* is the json object defined in the kernel file.
+    *info* is the json object defined in the kernel file.
 
     *form* is either "Iq" or "Iqxy".
@@ -206,32 +212,16 @@
     # If we are building the Iqxy kernel, we need to propagate qx,qy
     # parameters, otherwise we can
-    if form == "Iqxy":
-        qpars = IQXY_KERNEL
-    else:
-        qpars = IQ_KERNEL
-
+    dim = "2d" if is_2D else "1d"
+    fixed_pars = info['partype']['fixed-'+dim]
+    pd_pars = info['partype']['pd-'+dim]
+    vol_pars = info['partype']['volume']
+    q_pars = KERNEL_2D if is_2D else KERNEL_1D
+
+    # Build polydispersity loops
     depth = 4
     offset = ""
     loop_head = []
     loop_end = []
-    vol_pars = []
-    fixed_pars = []
-    pd_pars = []
-    fn_pars = []
-    for i,p in enumerate(meta['parameters']):
-        name = p[0]
-        ptype = p[4]
-        if ptype == "volume":
-            vol_pars.append(name)
-        elif ptype == "orientation":
-            if form != "Iqxy": continue  # no orientation for 1D kernels
-        elif ptype == "magnetic":
-            raise NotImplementedError("no magnetic parameters yet")
-        if name not in ['scale','background']: fn_pars.append(name)
-        if ptype == "":
-            fixed_pars.append(name)
-            continue
-        else:
-            pd_pars.append(name)
+    for name in pd_pars:
         subst = { 'name': name, 'offset': offset }
         loop_head.append(indent(LOOP_OPEN%subst, depth))
@@ -254,10 +244,16 @@
 
     # Define the inner loop function call
+    # The parameters to the f(q,p1,p2...) call should occur in the same
+    # order as given in the parameter info structure.  This may be different
+    # from the parameter order in the call to the kernel since the kernel
+    # call places all fixed parameters before all polydisperse parameters.
+    fq_pars = [p[0] for p in info['parameters'][len(COMMON_PARAMETERS):]
+               if p[0] in set(fixed_pars+pd_pars)]
     subst = {
         'weight_product': "*".join(p+"_w" for p in pd_pars),
         'volume_norm': volume_norm,
-        'fn': qpars['fn'],
-        'qcall': qpars['qcall'],
-        'pcall': ", ".join(fn_pars),
+        'fn': q_pars['fn'],
+        'qcall': q_pars['qcall'],
+        'pcall': ", ".join(fq_pars), # skip scale and background
         }
     loop_body = [indent(LOOP_BODY%subst, depth)]
@@ -280,7 +276,7 @@
     subst = {
         # kernel name is, e.g., cylinder_Iq
-        'name': "_".join((meta['name'], qpars['fn'])),
+        'name': kernel_name(info, is_2D),
         # to declare, e.g., global real q[],
-        'q_par_decl': qpars['q_par_decl'],
+        'q_par_decl': q_pars['q_par_decl'],
         # to declare, e.g., real sld, int Nradius, int Nlength
         'par_decl': par_decl,
@@ -288,5 +284,5 @@
         'pd_length': "+".join('N'+p for p in pd_pars),
         # the q initializers, e.g., real qi = q[i];
-        'qinit': qpars['qinit'],
+        'qinit': q_pars['qinit'],
         # the actual polydispersity loop
         'loops': loops,
@@ -295,6 +291,6 @@
     return kernel
 
-def make_partable(meta):
-    pars = meta['parameters']
+def make_partable(info):
+    pars = info['parameters']
     column_widths = [
         max(len(p[0]) for p in pars),
@@ -320,13 +316,13 @@
     return "\n".join(lines)
 
-def make_doc(kernelfile, meta, doc):
-    doc = doc%{'parameters': make_partable(meta)}
+def make_doc(kernelfile, info, doc):
+    doc = doc%{'parameters': make_partable(info)}
     return doc
 
-def make_model(kernelfile, meta, source):
-    kernel_Iq = make_kernel(meta, "Iq")
-    kernel_Iqxy = make_kernel(meta, "Iqxy")
+def make_model(kernelfile, info, source):
+    kernel_Iq = make_kernel(info, is_2D=False)
+    kernel_Iqxy = make_kernel(info, is_2D=True)
     path = os.path.dirname(kernelfile)
-    extra = [open("%s/%s"%(path,f)).read() for f in meta['include']]
+    extra = [open("%s/%s"%(path,f)).read() for f in info['include']]
     kernel = "\n\n".join([KERNEL_HEADER]+extra+[source, kernel_Iq, kernel_Iqxy])
     return kernel
@@ -339,20 +335,18 @@
     if len(parts) != 3:
         raise ValueError("PARAMETERS block missing from %r"%kernelfile)
-    meta_source = parts[1].strip()
+    info_source = parts[1].strip()
     try:
-        meta = relaxed_loads(meta_source)
+        info = relaxed_loads(info_source)
     except:
         print "in json text:"
         print "\n".join("%2d: %s"%(i+1,s)
-                        for i,s in enumerate(meta_source.split('\n')))
+                        for i,s in enumerate(info_source.split('\n')))
         raise
         #raise ValueError("PARAMETERS block could not be parsed from %r"%kernelfile)
-    meta['parameters'] = COMMON_PARAMETERS + meta['parameters']
-    meta['filename'] = kernelfile
 
     # select documentation out of the source file
     parts = source.split("DOCUMENTATION")
     if len(parts) == 3:
-        doc = make_doc(kernelfile, meta, parts[1].strip())
+        doc = make_doc(kernelfile, info, parts[1].strip())
     elif len(parts) == 1:
         raise ValueError("DOCUMENTATION block is missing from %r"%kernelfile)
@@ -360,5 +354,57 @@
         raise ValueError("DOCUMENTATION block incorrect from %r"%kernelfile)
 
-    return source, meta, doc
+    return source, info, doc
+
+def categorize_parameters(pars):
+    """
+    Build parameter categories out of the the parameter definitions.
+
+    Returns a dictionary of categories.
+
+    The function call sequence consists of q inputs and the return vector,
+    followed by the loop value/weight vector, followed by the values for
+    the non-polydisperse parameters, followed by the lengths of the
+    polydispersity loops.  To construct the call for 1D models, the
+    categories *fixed-1d* and *pd-1d* list the names of the parameters
+    of the non-polydisperse and the polydisperse parameters respectively.
+    Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models.
+    The *pd-rel* category is a set of those parameters which give
+    polydispersitiy as a portion of the value (so a 10% length dispersity
+    would use a polydispersity value of 0.1) rather than absolute
+    dispersity such as an angle plus or minus 15 degrees.
+
+    The *volume* category lists the volume parameters in order for calls
+    to volume within the kernel (used for volume normalization) and for
+    calls to ER and VR for effective radius and volume ratio respectively.
+
+    The *orientation* and *magnetic* categories list the orientation and
+    magnetic parameters.  These are used by the sasview interface.  The
+    blank category is for parameters such as scale which don't have any
+    other marking.
+    """
+    partype = {
+        'volume': [], 'orientation': [], 'magnetic': [], '': [],
+        'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [],
+        'pd-rel': set(),
+    }
+
+    for p in pars:
+        name,ptype = p[0],p[4]
+        if ptype == 'volume':
+            partype['pd-1d'].append(name)
+            partype['pd-2d'].append(name)
+            partype['pd-rel'].add(name)
+        elif ptype == 'magnetic':
+            partype['fixed-2d'].append(name)
+        elif ptype == 'orientation':
+            partype['pd-2d'].append(name)
+        elif ptype == '':
+            partype['fixed-1d'].append(name)
+            partype['fixed-2d'].append(name)
+        else:
+            raise ValueError("unknown parameter type %r"%ptype)
+        partype[ptype].append(name)
+
+    return partype
 
 def make(kernelfile):
@@ -370,59 +416,12 @@
     """
     #print kernelfile
-    source, meta, doc = parse_file(kernelfile)
-    doc = make_doc(kernelfile, meta, doc)
-    model = make_model(kernelfile, meta, source)
-    return model, meta, doc
-
-
-
-# Convert from python float to C float or double, depending on dtype
-FLOAT_CONVERTER = {
-    F32: np.float32,
-    F64: np.float64,
-    }
-
-def kernel_name(meta, is_2D):
-    return meta['name'] + "_" + ("Iqxy" if is_2D else "Iq")
-
-
-def kernel_pars(pars, par_info, is_2D, dtype=F32):
-    """
-    Convert parameter dictionary into arguments for the kernel.
-
-    *pars* is a dictionary of *{ name: value }*, with *name_pd* for the
-    polydispersity width, *name_pd_n* for the number of pd steps, and
-    *name_pd_nsigma* for the polydispersity limits.
-
-    *par_info* is the parameter info structure from the kernel metadata.
-
-    *is_2D* is True if the dataset represents 2D data, with the corresponding
-    orientation parameters.
-
-    *dtype* is F32 or F64, the numpy single and double precision floating
-    point dtypes.  These should not be the strings.
-    """
-    from .weights import GaussianDispersion
-    real = np.float32 if dtype == F32 else np.float64
-    fixed = []
-    parts = []
-    for p in par_info['parameters']:
-        name, ptype = p[0],p[4]
-        value = pars[name]
-        if ptype == "":
-            fixed.append(real(value))
-        elif is_2D or ptype != "orientation":
-            limits, width = p[3], pars[name+'_pd']
-            n, nsigma = pars[name+'_pd_n'], pars[name+'_pd_nsigma']
-            relative = (ptype != "orientation")
-            dist = GaussianDispersion(int(n), width, nsigma)
-            # Make sure that weights are normalized to peaks at 1 so that
-            # the tolerance term can be used properly on truncated distributions
-            v,w = dist.get_weights(value, limits[0], limits[1], relative)
-            parts.append((v, w/w.max()))
-    loops = np.hstack(parts)
-    loops = np.ascontiguousarray(loops.T, dtype).flatten()
-    loopsN = [np.uint32(len(p[0])) for p in parts]
-    return fixed, loops, loopsN
+    source, info, doc = parse_file(kernelfile)
+    info['filename'] = kernelfile
+    info['parameters'] = COMMON_PARAMETERS + info['parameters']
+    info['partype'] = categorize_parameters(info['parameters'])
+    info['limits'] = dict((p[0],p[3]) for p in info['parameters'])
+    doc = make_doc(kernelfile, info, doc)
+    model = make_model(kernelfile, info, source)
+    return model, info, doc
 
 
@@ -437,5 +436,5 @@
 def demo():
     from os.path import join as joinpath, dirname
-    c, meta, doc = make_model(joinpath(dirname(__file__), "models", "cylinder.c"))
+    c, info, doc = make_model(joinpath(dirname(__file__), "models", "cylinder.c"))
     #print doc
     #print c
Index: sasmodels/gpu.py
===================================================================
--- sasmodels/gpu.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/gpu.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -45,4 +45,12 @@
 #define real double
 """
+
+# The max loops number is limited by the amount of local memory available
+# on the device.  You don't want to make this value too big because it will
+# waste resources, nor too small because it may interfere with users trying
+# to do their polydispersity calculations.  A value of 1024 should be much
+# larger than necessary given that cost grows as npts^k where k is the number
+# of polydisperse parameters.
+MAX_LOOPS = 1024
 
 ENV = None
@@ -96,6 +104,5 @@
     dtype = np.dtype(dtype)
     if dtype==F64 and not all(has_double(d) for d in context.devices):
-        warnings.warn(RuntimeWarning("Double precision not support; using single precision instead"))
-        dtype = F32
+        raise RuntimeError("Double precision not supported for devices")
 
     header = F64_DEFS if dtype == F64 else F32_DEFS
@@ -104,5 +111,5 @@
         header += "#define USE_SINCOS\n"
     program  = cl.Program(context, header+source).build()
-    return program, dtype
+    return program
 
 
@@ -125,4 +132,17 @@
         self.boundary = max(d.min_data_type_align_size
                             for d in self.context.devices)
+        self.has_double = all(has_double(d) for d in self.context.devices)
+        self.compiled = {}
+
+    def compile_program(self, name, source, dtype):
+        if name not in self.compiled:
+            #print "compiling",name
+            self.compiled[name] = compile_model(self.context, source, dtype)
+        return self.compiled[name]
+
+    def release_program(self, name):
+        if name in self.compiled:
+            self.compiled[name].release()
+            del self.compiled[name]
 
 class GpuModel(object):
@@ -130,5 +150,5 @@
     GPU wrapper for a single model.
 
-    *source* and *meta* are the model source and interface as returned
+    *source* and *info* are the model source and interface as returned
     from :func:`gen.make`.
 
@@ -138,19 +158,31 @@
     is an optional extension which may not be available on all devices.
     """
-    def __init__(self, source, meta, dtype=F32):
-        context = environment().context
-        self.meta = meta
-        self.program, self.dtype = compile_model(context, source, dtype)
-
-        #TODO: need better interface
-        self.PARS = dict((p[0],p[2]) for p in meta['parameters'])
-        self.PD_PARS = [p[0] for p in meta['parameters'] if p[4] != ""]
-
-    def __call__(self, input, cutoff=1e-5):
+    def __init__(self, source, info, dtype=F32):
+        self.info = info
+        self.source = source
+        self.dtype = dtype
+        self.program = None # delay program creation
+
+    def __getstate__(self):
+        state = self.__dict__.copy()
+        state['program'] = None
+        return state
+
+    def __setstate__(self, state):
+        self.__dict__ = state.copy()
+
+    def __call__(self, input):
         if self.dtype != input.dtype:
             raise TypeError("data and kernel have different types")
-        kernel_name = gen.kernel_name(self.meta, input.is_2D)
+        if self.program is None:
+            self.program = environment().compile_program(self.info['name'],self.source, self.dtype)
+        kernel_name = gen.kernel_name(self.info, input.is_2D)
         kernel = getattr(self.program, kernel_name)
-        return GpuKernel(kernel, self.meta, input, cutoff)
+        return GpuKernel(kernel, self.info, input)
+
+    def release(self):
+        if self.program is not None:
+            environment().release_program(self.info['name'])
+            self.program = None
 
     def make_input(self, q_vectors):
@@ -210,15 +242,18 @@
 
 class GpuKernel(object):
-    def __init__(self, kernel, meta, input, cutoff):
-        env = environment()
-
-        self.cutoff = cutoff
+    def __init__(self, kernel, info, input):
         self.input = input
         self.kernel = kernel
-        self.meta = meta
+        self.info = info
+        self.res = np.empty(input.nq, input.dtype)
+        dim = '2d' if input.is_2D else '1d'
+        self.fixed_pars = info['partype']['fixed-'+dim]
+        self.pd_pars = info['partype']['pd-'+dim]
 
         # Inputs and outputs for each kernel call
+        # Note: res may be shorter than res_b if global_size != nq
+        env = environment()
         self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                  1024*input.dtype.itemsize)
+                                  MAX_LOOPS*input.dtype.itemsize)
                         for _ in env.queues]
         self.res_b = [cl.Buffer(env.context, mf.READ_WRITE,
@@ -226,28 +261,26 @@
                       for _ in env.queues]
 
-        # Note: may be shorter than res_b if global_size != nq
-        self.res = np.empty(input.nq, input.dtype)
-
-        # Determine the set of fixed and polydisperse parameters
-        self.fixed_pars = [p[0] for p in meta['parameters'] if p[4] == ""]
-        self.pd_pars = [p for p in meta['parameters']
-               if p[4]=="volume" or (p[4]=="orientation" and input.is_2D)]
-
-    def eval(self, pars):
+
+    def __call__(self, pars, pd_pars, cutoff=1e-5):
+        real = np.float32 if self.input.dtype == F32 else np.float64
+        fixed = [real(p) for p in pars]
+        cutoff = real(cutoff)
+        loops = np.hstack(pd_pars)
+        loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten()
+        loops_N = [np.uint32(len(p[0])) for p in pd_pars]
+
+        #print "opencl eval",pars
+        if len(loops) > MAX_LOOPS:
+            raise ValueError("too many polydispersity points")
         device_num = 0
         res_bi = self.res_b[device_num]
         queuei = environment().queues[device_num]
-        fixed, loops, loop_n = \
-            gen.kernel_pars(pars, self.meta, self.input.is_2D, dtype=self.input.dtype)
+        loops_bi = self.loops_b[device_num]
         loops_l = cl.LocalMemory(len(loops.data))
-        real = np.float32 if self.input.dtype == F32 else np.float64
-        cutoff = real(self.cutoff)
-
-        loops_bi = self.loops_b[device_num]
         cl.enqueue_copy(queuei, loops_bi, loops)
         #ctx = environment().context
         #loops_bi = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops)
-        pars = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + loop_n
-        self.kernel(queuei, self.input.global_size, None, *pars)
+        args = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + loops_N
+        self.kernel(queuei, self.input.global_size, None, *args)
         cl.enqueue_copy(queuei, self.res, res_bi)
 
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/models/cylinder.c	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -140,14 +140,13 @@
 DOCUMENTATION END
 */
-real volume(real radius, real length);
+real form_volume(real radius, real length);
 real Iq(real q, real sld, real solvent_sld, real radius, real length);
 real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi);
 
 
-real volume(real radius, real length)
+real form_volume(real radius, real length)
 {
     return M_PI*radius*radius*length;
 }
-
 
 real Iq(real q,
@@ -170,5 +169,5 @@
     // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
     // The additional volume factor is for polydisperse volume normalization.
-    const real s = (sld - solvent_sld) * volume(radius, length);
+    const real s = (sld - solvent_sld) * form_volume(radius, length);
     return REAL(1.0e-4) * form * s * s;
 }
@@ -203,10 +202,10 @@
     const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg);
     const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg);
-    const real form = REAL(4.0)*bj*bj*si*si;
+    const real form = bj*bj*si*si;
 
     // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
     // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
     // The additional volume factor is for polydisperse volume normalization.
-    const real s = (sld - solvent_sld) * volume(radius, length);
+    const real s = (sld - solvent_sld) * form_volume(radius, length);
     return REAL(1.0e-4) * form * s * s; // * correction;
 }
Index: sasmodels/models/cylinder_clone.c
===================================================================
--- sasmodels/models/cylinder_clone.c	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
+++ sasmodels/models/cylinder_clone.c	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -0,0 +1,212 @@
+/* PARAMETERS
+{
+name: "CylinderModel",
+title: "Cylinder with uniform scattering length density",
+include: [ "lib/J1.c", "lib/gauss76.c", "lib/cylkernel.c"],
+parameters: [
+   // [ "name", "units", default, [lower, upper], "type", "description" ],
+   [ "sldCyl", "1e-6/Ang^2", 4e-6, [-Infinity,Infinity], "",
+     "Cylinder scattering length density" ],
+   [ "sldSolv", "1e-6/Ang^2", 1e-6, [-Infinity,Infinity], "",
+     "Solvent scattering length density" ],
+   [ "radius", "Ang",  20, [0, Infinity], "volume",
+     "Cylinder radius" ],
+   [ "length", "Ang",  400, [0, Infinity], "volume",
+     "Cylinder length" ],
+   [ "cyl_theta", "degrees", 60, [-Infinity, Infinity], "orientation",
+     "In plane angle" ],
+   [ "cyl_phi", "degrees", 60, [-Infinity, Infinity], "orientation",
+     "Out of plane angle" ],
+],
+description: "f(q)= 2*(sldCyl - sldSolv)*V*sin(qLcos(alpha/2))/[qLcos(alpha/2)]*J1(qRsin(alpha/2))/[qRsin(alpha)]",
+}
+PARAMETERS END
+
+DOCUMENTATION
+.. _CylinderModel:
+
+CylinderModel
+=============
+
+This model provides the form factor for a right circular cylinder with uniform
+scattering length density. The form factor is normalized by the particle volume.
+
+For information about polarised and magnetic scattering, click here_.
+
+Definition
+----------
+
+The output of the 2D scattering intensity function for oriented cylinders is
+given by (Guinier, 1955)
+
+.. math::
+
+    P(q,\alpha) = \frac{\text{scale}}{V}f^2(q) + \text{bkg}
+
+where
+
+.. math::
+
+    f(q) = 2 (\Delta \rho) V
+           \frac{\sin (q L/2 \cos \alpha)}{q L/2 \cos \alpha}
+           \frac{J_1 (q r \sin \alpha)}{q r \sin \alpha}
+
+and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V$
+is the volume of the cylinder, $L$ is the length of the cylinder, $r$ is the
+radius of the cylinder, and $d\rho$ (contrast) is the scattering length density
+difference between the scatterer and the solvent. $J_1$ is the first order
+Bessel function.
+
+To provide easy access to the orientation of the cylinder, we define the
+axis of the cylinder using two angles $\theta$ and $\phi$. Those angles
+are defined in Figure :num:`figure #CylinderModel-orientation`.
+
+.. _CylinderModel-orientation:
+
+.. figure:: img/image061.JPG
+
+    Definition of the angles for oriented cylinders.
+
+.. figure:: img/image062.JPG
+
+    Examples of the angles for oriented pp against the detector plane.
+
+NB: The 2nd virial coefficient of the cylinder is calculated based on the
+radius and length values, and used as the effective radius for $S(Q)$
+when $P(Q) \cdot S(Q)$ is applied.
+
+The returned value is scaled to units of |cm^-1| and the parameters of
+the CylinderModel are the following:
+
+%(parameters)s
+
+The output of the 1D scattering intensity function for randomly oriented
+cylinders is then given by
+
+.. math::
+
+    P(q) = \frac{\text{scale}}{V}
+        \int_0^{\pi/2} f^2(q,\alpha) \sin \alpha d\alpha + \text{background}
+
+The *theta* and *phi* parameters are not used for the 1D output. Our
+implementation of the scattering kernel and the 1D scattering intensity
+use the c-library from NIST.
+
+Validation of the CylinderModel
+-------------------------------
+
+Validation of our code was done by comparing the output of the 1D model
+to the output of the software provided by the NIST (Kline, 2006).
+Figure :num:`figure #CylinderModel-compare` shows a comparison of
+the 1D output of our model and the output of the NIST software.
+
+.. _CylinderModel-compare:
+
+.. figure:: img/image065.JPG
+
+    Comparison of the SasView scattering intensity for a cylinder with the
+    output of the NIST SANS analysis software.
+    The parameters were set to: *Scale* = 1.0, *Radius* = 20 |Ang|,
+    *Length* = 400 |Ang|, *Contrast* = 3e-6 |Ang^-2|, and
+    *Background* = 0.01 |cm^-1|.
+
+In general, averaging over a distribution of orientations is done by
+evaluating the following
+
+.. math::
+
+    P(q) = \int_0^{\pi/2} d\phi
+        \int_0^\pi p(\theta, \phi) P_0(q,\alpha) \sin \theta d\theta
+
+
+where $p(\theta,\phi)$ is the probability distribution for the orientation
+and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented
+system. Since we have no other software to compare the implementation of
+the intensity for fully oriented cylinders, we can compare the result of
+averaging our 2D output using a uniform distribution $p(\theta, \phi) = 1.0$.
+Figure :num:`figure #CylinderModel-crosscheck` shows the result of
+such a cross-check.
+
+.. _CylinderModel-crosscheck:
+
+.. figure:: img/image066.JPG
+
+    Comparison of the intensity for uniformly distributed cylinders
+    calculated from our 2D model and the intensity from the NIST SANS
+    analysis software.
+    The parameters used were: *Scale* = 1.0, *Radius* = 20 |Ang|,
+    *Length* = 400 |Ang|, *Contrast* = 3e-6 |Ang^-2|, and
+    *Background* = 0.0 |cm^-1|.
+
+DOCUMENTATION END
+*/
+real form_volume(real radius, real length);
+real Iq(real q, real sld, real solvent_sld, real radius, real length);
+real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi);
+
+
+real form_volume(real radius, real length)
+{
+    return M_PI*radius*radius*length;
+}
+
+real Iq(real q,
+    real sldCyl,
+    real sldSolv,
+    real radius,
+    real length)
+{
+    const real h = REAL(0.5)*length;
+    real summ = REAL(0.0);
+    for (int i=0; i<76 ;i++) {
+        //const real zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+        const real zi = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);
+        summ += Gauss76Wt[i] * CylKernel(q, radius, h, zi);
+    }
+    //const real form = (uplim-lolim)/2.0*summ;
+    const real form = summ * M_PI_4;
+
+    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
+    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
+    // The additional volume factor is for polydisperse volume normalization.
+    const real s = (sldCyl - sldSolv) * form_volume(radius, length);
+    return REAL(1.0e8) * form * s * s;
+}
+
+
+real Iqxy(real qx, real qy,
+    real sldCyl,
+    real sldSolv,
+    real radius,
+    real length,
+    real cyl_theta,
+    real cyl_phi)
+{
+    real sn, cn; // slots to hold sincos function output
+
+    // Compute angle alpha between q and the cylinder axis
+    SINCOS(cyl_theta*M_PI_180, sn, cn);
+    // # The following correction factor exists in sasview, but it can't be
+    // # right, so we are leaving it out for now.
+    // const real correction = fabs(cn)*M_PI_2;
+    const real q = sqrt(qx*qx+qy*qy);
+    const real cos_val = cn*cos(cyl_phi*M_PI_180)*(qx/q) + sn*(qy/q);
+    const real alpha = acos(cos_val);
+
+    // The following is CylKernel() / sin(alpha), but we are doing it in place
+    // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit
+    // faster since we don't mulitply and divide sin(alpha).
+    SINCOS(alpha, sn, cn);
+    const real besarg = q*radius*sn;
+    const real siarg = REAL(0.5)*q*length*cn;
+    // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
+    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg);
+    const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg);
+    const real form = bj*bj*si*si;
+
+    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
+    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
+    // The additional volume factor is for polydisperse volume normalization.
+    const real s = (sldCyl - sldSolv) * form_volume(radius, length);
+    return REAL(1.0e8) * form * s * s; // * correction;
+}
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
+++ sasmodels/models/ellipsoid.c	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -0,0 +1,85 @@
+/* PARAMETERS
+{
+name: "ellipsoid",
+title: "Ellipsoid with uniform scattering length density",
+include: [ "lib/gauss76.c" ],
+parameters: [
+   // [ "name", "units", default, [lower, upper], "type", "description" ],
+   [ "sld", "1e-6/Ang^2", 4, [-Infinity,Infinity], "",
+     "Cylinder scattering length density" ],
+   [ "solvent_sld", "1e-6/Ang^2", 1, [-Infinity,Infinity], "",
+     "Solvent scattering length density" ],
+   [ "a", "Ang",  20, [0, Infinity], "volume",
+     "Cylinder radius" ],
+   [ "b", "Ang",  20, [0, Infinity], "volume",
+     "Cylinder length" ],
+   [ "theta", "degrees", 60, [-Infinity, Infinity], "orientation",
+     "In plane angle" ],
+   [ "phi", "degrees", 60, [-Infinity, Infinity], "orientation",
+     "Out of plane angle" ],
+],
+}
+PARAMETERS END
+
+DOCUMENTATION
+.. _EllipseModel:
+
+DOCUMENTATION END
+*/
+
+real form_volume(real a, real b);
+real Iq(real qx, real qy, real sld, real solvent_sld, real a, real b);
+real Iqxy(real qx, real qy, real sld, real solvent_sld, real a, real b, real theta, real phi);
+
+real form_volume(real a, real b)
+{
+    return REAL(1.333333333333333)*M_PI_2*a*b*b;
+}
+
+real ellipsoid_kernel(double q, double b, double a, double dum)
+{
+    real sn, cn;
+    const real nu = a/b;
+    const real arg = q * b * sqrt(REAL(1.0)+(dum*dum*(nu*nu--REAL(1.0))));
+    SINCOS(arg, sn, cn);
+    const real f = (arg==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-arg*cn)/(arg*arg*arg);
+    return f*f;
+}
+
+real Iq(real q,
+    real sld,
+    real solvent_sld,
+    real a,
+    real b)
+{
+    real summ = REAL(0.0);
+    for (int i=0;i<76;i++) {
+        //const real zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+        zi = ( Gauss76Z[i] + REAL(1.0))/REAL(2.0);
+        summ += Gauss76Wt[i] * ellipsoid_kernel(q, b, a, zi);
+    }
+    //const real form = (uplim-lolim)/2.0*summ;
+    const real form = REAL(0.5)*summ
+    const real s = (sld - sld_solvent) * form_volume(a, b);
+    return REAL(1.0e-4) * form * s * s;
+}
+
+real Iqxy(real qx, real qy,
+    real sld,
+    real solvent_sld,
+    real a,
+    real b,
+    real theta,
+    real phi)
+{
+    real sn, cn;
+
+    const real q = sqrt(qx*qx + qy*qy);
+    SINCOS(theta*M_PI_180, sn, cn);
+    const real cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
+    const real form = ellipsoid_kernel(q, b, a, cos_val);
+    const real s = (sld - solvent_sld) * form_volume(a, b);
+
+    return REAL(1.0e-4) * form * s * s;
+}
+
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
+++ sasmodels/sasview_model.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -0,0 +1,325 @@
+import math
+from copy import deepcopy
+
+import numpy as np
+
+def make_class(name, class_name=None):
+    from .core import opencl_model
+    if class_name is None:
+        class_name = "".join(part.capitalize() for part in name.split('_')+['model'])
+    model =  opencl_model(name)
+    class ConstructedModel(SasviewModel):
+        kernel = model
+        def __init__(self, multfactor=1):
+            if self.kernel is None:
+                self.__class__.kernel = opencl_model(name)
+            SasviewModel.__init__(self, self.kernel)
+    ConstructedModel.__name__ = class_name
+    return ConstructedModel
+
+class SasviewModel(object):
+    """
+    Sasview wrapper for opencl/ctypes model.
+    """
+    def __init__(self, model):
+        """Initialization"""
+        self._model = model
+
+        self.name = model.info['name']
+        self.description = model.info['description']
+        self.category = None
+        self.multiplicity_info = None
+        self.is_multifunc = False
+
+        ## interpret the parameters
+        ## TODO: reorganize parameter handling
+        self.details = dict()
+        self.params = dict()
+        self.dispersion = dict()
+        partype = model.info['partype']
+        for name,units,default,limits,ptype,description in model.info['parameters']:
+            self.params[name] = default
+            self.details[name] = [units]+limits
+
+        for name in partype['pd-2d']:
+            self.dispersion[name] = {
+                'width': 0,
+                'npts': 35,
+                'nsigma': 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']]
+        self.non_fittable = []
+
+        ## independent parameter name and unit [string]
+        self.input_name = model.info.get("input_name","Q")
+        self.input_unit = model.info.get("input_unit","A^{-1}")
+        self.output_name = model.info.get("output_name","Intensity")
+        self.output_unit = model.info.get("output_unit","cm^{-1}")
+
+        ## _persistency_dict is used by sans.perspectives.fitting.basepage
+        ## to store dispersity reference.
+        ## TODO: _persistency_dict to persistency_dict throughout sasview
+        self._persistency_dict = {}
+
+        ## New fields introduced for opencl rewrite
+        self.cutoff = 1e-5
+
+    def __str__(self):
+        """
+        :return: string representation
+        """
+        return self.name
+
+    def is_fittable(self, par_name):
+        """
+        Check if a given parameter is fittable or not
+
+        :param par_name: the parameter name to check
+        """
+        return par_name.lower() in self.fixed
+        #For the future
+        #return self.params[str(par_name)].is_fittable()
+
+
+    def getProfile(self):
+        """
+        Get SLD profile
+
+        : return: (z, beta) where z is a list of depth of the transition points
+                beta is a list of the corresponding SLD values
+        """
+        return None, None
+
+    def setParam(self, name, value):
+        """
+        Set the value of a model parameter
+
+        :param name: name of the parameter
+        :param value: value of the parameter
+
+        """
+        # Look for dispersion parameters
+        toks = name.split('.')
+        if len(toks)==2:
+            for item in self.dispersion.keys():
+                if item.lower()==toks[0].lower():
+                    for par in self.dispersion[item]:
+                        if par.lower() == toks[1].lower():
+                            self.dispersion[item][par] = value
+                            return
+        else:
+            # Look for standard parameter
+            for item in self.params.keys():
+                if item.lower()==name.lower():
+                    self.params[item] = value
+                    return
+
+        raise ValueError, "Model does not contain parameter %s" % name
+
+    def getParam(self, name):
+        """
+        Set the value of a model parameter
+
+        :param name: name of the parameter
+
+        """
+        # Look for dispersion parameters
+        toks = name.split('.')
+        if len(toks)==2:
+            for item in self.dispersion.keys():
+                if item.lower()==toks[0].lower():
+                    for par in self.dispersion[item]:
+                        if par.lower() == toks[1].lower():
+                            return self.dispersion[item][par]
+        else:
+            # Look for standard parameter
+            for item in self.params.keys():
+                if item.lower()==name.lower():
+                    return self.params[item]
+
+        raise ValueError, "Model does not contain parameter %s" % name
+
+    def getParamList(self):
+        """
+        Return a list of all available parameters for the model
+        """
+        list = self.params.keys()
+        # WARNING: Extending the list with the dispersion parameters
+        list.extend(self.getDispParamList())
+        return list
+
+    def getDispParamList(self):
+        """
+        Return a list of all available parameters for the model
+        """
+        list = []
+
+        for item in self.dispersion.keys():
+            for p in self.dispersion[item].keys():
+                if p not in ['type']:
+                    list.append('%s.%s' % (item.lower(), p.lower()))
+
+        return list
+
+    def clone(self):
+        """ Return a identical copy of self """
+        return deepcopy(self)
+
+    def run(self, x=0.0):
+        """
+        Evaluate the model
+
+        :param x: input q, or [q,phi]
+
+        :return: scattering function P(q)
+
+        **DEPRECATED**: use calculate_Iq instead
+        """
+        if isinstance(x, (list,tuple)):
+            q, phi = x
+            return self.calculate_Iq([q * math.cos(phi)],
+                                     [q * math.sin(phi)])[0]
+        else:
+            return self.calculate_Iq([float(x)])[0]
+
+
+    def runXY(self, x=0.0):
+        """
+        Evaluate the model in cartesian coordinates
+
+        :param x: input q, or [qx, qy]
+
+        :return: scattering function P(q)
+
+        **DEPRECATED**: use calculate_Iq instead
+        """
+        if isinstance(x, (list,tuple)):
+            return self.calculate_Iq([float(x[0])],[float(x[1])])[0]
+        else:
+            return self.calculate_Iq([float(x)])[0]
+
+    def evalDistribution(self, qdist):
+        """
+        Evaluate a distribution of q-values.
+
+        * For 1D, a numpy array is expected as input: ::
+
+            evalDistribution(q)
+
+          where q is a numpy array.
+
+        * For 2D, a list of numpy arrays are expected: [qx,qy],
+          with 1D arrays::
+
+              qx = [ qx[0], qx[1], qx[2], ....]
+
+          and::
+
+              qy = [ qy[0], qy[1], qy[2], ....]
+
+        Then get ::
+
+            q = numpy.sqrt(qx^2+qy^2)
+
+        that is a qr in 1D array::
+
+            q = [q[0], q[1], q[2], ....]
+
+
+        :param qdist: ndarray of scalar q-values or list [qx,qy] where qx,qy are 1D ndarrays
+        """
+        if isinstance(qdist, (list,tuple)):
+            # Check whether we have a list of ndarrays [qx,qy]
+            qx, qy = qdist
+            return self.calculate_Iq(qx, qy)
+
+        elif isinstance(qdist, np.ndarray):
+            # We have a simple 1D distribution of q-values
+            return self.calculate_Iq(qdist)
+
+        else:
+            raise TypeError("evalDistribution expects q or [qx, qy], not %r"%type(qdist))
+
+    def calculate_Iq(self, *args):
+        q_vectors = [np.asarray(q) for q in args]
+        fn = self._model(self._model.make_input(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.input.release()
+        fn.release()
+        return result
+
+    def calculate_ER(self):
+        """
+        Calculate the effective radius for P(q)*S(q)
+
+        :return: the value of the effective radius
+        """
+        ER = self._model.info.get('ER', None)
+        if ER is None:
+            return 1.0
+        else:
+            vol_pars = self._model.info['partype']['volume']
+            values, weights = self._dispersion_mesh(vol_pars)
+            fv = ER(*values)
+            return np.sum(weights*fv) / np.sum(weights)
+
+    def calculate_VR(self):
+        """
+        Calculate the volf ratio for P(q)*S(q)
+
+        :return: the value of the volf ratio
+        """
+        VR = self._model.info.get('ER', None)
+        if VR is None:
+            return 1.0
+        else:
+            vol_pars = self._model.info['partype']['volume']
+            values, weights = self._dispersion_mesh(vol_pars)
+            whole,part = VR(*values)
+            return np.sum(weights*part)/np.sum(weights*whole)
+
+    def set_dispersion(self, parameter, dispersion):
+        """
+        Set the dispersion object for a model parameter
+
+        :param parameter: name of the parameter [string]
+        :param dispersion: dispersion object of type Dispersion
+        """
+        if parameter.lower() in (s.name.lower() for s in self.params.keys()):
+            self.dispersion[parameter] = dispersion.get_pars()
+        else:
+            raise ValueError("%r is not a dispersity or orientation parameter")
+
+    def _dispersion_mesh(self, pars):
+        """
+        Create a mesh grid of dispersion parameters and weights.
+
+        Returns [p1,p2,...],w where pj is a vector of values for parameter j
+        and w is a vector containing the products for weights for each
+        parameter set in the vector.
+        """
+        values, weights = zip(*[self._get_weights(p) for p in pars])
+        values = [v.flatten() for v in np.meshgrid(*values)]
+        weights = np.vstack([v.flatten() for v in np.meshgrid(*weights)])
+        weights = np.prod(weights, axis=1)
+        return values, weights
+
+    def _get_weights(self, par):
+        from . import weights
+
+        relative = self._model.info['partype']['pd-rel']
+        limits = self._model.info['limits']
+        dis = self.dispersion[par]
+        v,w = weights.get_weights(
+            dis['type'], dis['npts'], dis['width'], dis['nsigma'],
+            self.params[par], limits[par], par in relative)
+        return v,w/w.max()
Index: sasmodels/weights.py
===================================================================
--- sasmodels/weights.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/weights.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
@@ -1,26 +1,128 @@
+# TODO: include dispersion docs with the disperser models
+from math import sqrt
 import numpy as np
+from scipy.special import gammaln
 
-class GaussianDispersion(object):
-    def __init__(self, npts=35, width=0, nsigmas=3): #number want, percent deviation, #standard deviations from mean
-        self.type = 'gaussian'
-        self.npts = npts
-        self.width = width
-        self.nsigmas = nsigmas
+class Dispersion(object):
+    """
+    Base dispersion object.
+
+    Subclasses should define *_weights(center, sigma, lb, ub)*
+    which returns the x points and their corresponding weights.
+    """
+    type = "base disperser"
+    default = dict(npts=35, width=0, nsigmas=3)
+    def __init__(self, npts=None, width=None, nsigmas=None):
+        self.npts = self.default['npts'] if npts is None else npts
+        self.width = self.default['width'] if width is None else width
+        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas
 
     def get_pars(self):
-        return self.__dict__
+        pars = {'type': self.type}
+        pars.update(self.__dict__)
+        return pars
 
-    def get_weights(self, center, min, max, relative):
-        """ *center* is the center of the distribution
-        *min*,*max* are the min, max allowed values
-        *relative* is True if the width is relative to the center instead of absolute
-        For polydispersity use relative.  For orientation parameters use absolute."""
+    def set_weights(self, values, weights):
+        raise RuntimeError("set_weights is only available for ArrayDispersion")
+
+    def get_weights(self, center, lb, ub, relative):
+        """
+        Return the weights for the distribution.
+
+        *center* is the center of the distribution
+
+        *lb*,*ub* are the min and max allowed values
+
+        *relative* is True if the distribution width is proportional to the
+        center value instead of absolute.  For polydispersity use relative.
+        For orientation parameters use absolute.
+        """
         npts, width, nsigmas = self.npts, self.width, self.nsigmas
-        sigma = width * center if relative else width
+        sigma = self.width * center if relative else self.width
         if sigma == 0:
-            return np.array([center],'d'), np.array([1.], 'd')
-        x = center + np.linspace(-nsigmas * sigma, +nsigmas * sigma, npts)
-        x = x[(x >= min) & (x <= max)]
+            return np.array([center], 'd'), np.array([1.], 'd')
+        return self._weights(center, sigma, lb, ub)
+
+    def _weights(self, center, sigma, lb, ub):
+        """actual work of computing the weights"""
+        raise NotImplementedError
+
+    def _linspace(self, center, sigma, lb, ub):
+        """helper function to provide linear spaced weight points within range"""
+        npts, nsigmas = self.npts, self.nsigmas
+        x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts)
+        x = x[(x >= lb) & (x <= ub)]
+        return x
+
+
+class GaussianDispersion(Dispersion):
+    type = "gaussian"
+    default = dict(npts=35, width=0, nsigmas=3)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma, lb, ub)
         px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
         return x, px
 
+
+class RectangleDispersion(Dispersion):
+    type = "rectangle"
+    default = dict(npts=35, width=0, nsigmas=1.70325)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma*sqrt(3.0), lb, ub)
+        px = np.ones_like(x)
+        return x, px
+
+
+class LogNormalDispersion(Dispersion):
+    type = "lognormal"
+    default = dict(npts=80, width=0, nsigmas=8)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
+        px = np.exp(-0.5*(np.log(x)-center)**2)/sigma**2/(x*sigma)
+        return x, px
+
+
+class SchulzDispersion(Dispersion):
+    type = "schulz"
+    default = dict(npts=80, width=0, nsigmas=8)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
+        R= x/center
+        z = (center/sigma)**2
+        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
+        px = np.exp(arg)
+        return x, px
+
+
+class ArrayDispersion(Dispersion):
+    type = "array"
+    default = dict(npts=35, width=0, nsigmas=1)
+    def __init__(self, npts=None, width=None, nsigmas=None):
+        Dispersion.__init__(self, npts, width, nsigmas)
+        self.values = np.array([0.], 'd')
+        self.weights = np.array([1.], 'd')
+
+    def set_weights(self, values, weights):
+        self.values = np.ascontiguousarray(values, 'd')
+        self.weights = np.ascontiguousarray(weights, 'd')
+        self.npts = len(values)
+
+    def _weights(self, center, sigma, lb, ub):
+        # TODO: interpolate into the array dispersion using npts
+        x = center + self.values*sigma
+        idx = (x>=lb)&(x<=ub)
+        x = x[idx]
+        px = self.weights[idx]
+        return x, px
+
+
+models = dict((d.type,d) for d in (
+    GaussianDispersion, RectangleDispersion,
+    ArrayDispersion, SchulzDispersion, LogNormalDispersion
+))
+
+def get_weights(disperser, n, width, nsigma, value, limits, relative):
+    cls = models[disperser]
+    obj = cls(n, width, nsigma)
+    v,w = obj.get_weights(value, limits[0], limits[1], relative)
+    return v,w
