Index: .gitignore
===================================================================
--- .gitignore	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ .gitignore	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,3 @@
+*.pyc
+*.cl
+*.so
Index: compare-new.py
===================================================================
--- compare-new.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ compare-new.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,292 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+
+import numpy as np
+
+from sasmodels.core import BumpsModel, fake_data2D, set_beam_stop, plot_data, \
+    tic, opencl_model, dll_model
+
+def sasview_model(modelname, **pars):
+    """
+    Load a sasview model given the model name.
+    """
+    modelname = modelname+"Model"
+    sans = __import__('sans.models.'+modelname)
+    ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None)
+    if ModelClass is None:
+        raise ValueError("could not find model %r in sans.models"%modelname)
+    model = ModelClass()
+
+    for k,v in pars.items():
+        if k.endswith("_pd"):
+            model.dispersion[k[:-3]]['width'] = v
+        elif k.endswith("_pd_n"):
+            model.dispersion[k[:-5]]['npts'] = v
+        elif k.endswith("_pd_nsigma"):
+            model.dispersion[k[:-10]]['nsigmas'] = v
+        else:
+            model.setParam(k, v)
+    return model
+
+
+def compare(Ncpu, cpuname, cpupars, Ngpu, gpuname, gpupars):
+
+    #from sasmodels.core import load_data
+    #data = load_data('December/DEC07098.DAT')
+    data = fake_data2D(np.linspace(-0.05, 0.05, 128))
+    set_beam_stop(data, 0.004)
+
+    if Ngpu > 0:
+        gpumodel = opencl_model(gpuname, dtype='single')
+        model = BumpsModel(data, gpumodel, **gpupars)
+        toc = tic()
+        for i in range(Ngpu):
+            #pars['scale'] = np.random.rand()
+            model.update()
+            gpu = model.theory()
+        gpu_time = toc()*1000./Ngpu
+        print "ocl t=%.1f ms"%gpu_time
+        #print max(gpu), min(gpu)
+
+    if 0 and Ncpu > 0: # Hack to compare ctypes vs. opencl
+        dllmodel = dll_model(gpuname)
+        model = BumpsModel(data, dllmodel, **gpupars)
+        toc = tic()
+        for i in range(Ncpu):
+            model.update()
+            cpu = model.theory()
+        cpu_time = toc()*1000./Ncpu
+        print "dll t=%.1f ms"%cpu_time
+
+    elif 1: # Hack to check new vs old for GpuCylinder
+        from Models.code_cylinder_f import GpuCylinder as oldgpu
+        from sasmodel import SasModel
+        oldmodel = SasModel(data, oldgpu, dtype='single', **cpupars)
+        toc = tic()
+        for i in range(Ngpu):
+            oldmodel.update()
+            cpu = oldmodel.theory()
+        cpu_time = toc()*1000./Ngpu
+        print "old t=%.1f ms"%cpu_time
+
+    elif Ncpu > 0:
+        cpumodel = sasview_model(cpuname, **cpupars)
+        toc = tic()
+        for i in range(Ncpu):
+            cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data])
+        cpu_time = toc()*1000./Ncpu
+        print "sasview t=%.1f ms"%cpu_time
+
+    if Ngpu > 0 and Ncpu > 0:
+        print "gpu/cpu", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu))
+        #cpu *= max(gpu/cpu)
+        abserr = (gpu - cpu)
+        relerr = (gpu - cpu)/cpu
+        print "max(|ocl-omp|)", max(abs(abserr[model.index]))
+        print "max(|(ocl-omp)/ocl|)", max(abs(relerr[model.index]))
+    #return
+
+    import matplotlib.pyplot as plt
+    if Ncpu > 0:
+        if Ngpu > 0: plt.subplot(131)
+        plot_data(data, cpu)
+        plt.title("omp t=%.1f ms"%cpu_time)
+    if Ngpu > 0:
+        if Ncpu > 0: plt.subplot(132)
+        plot_data(data, gpu)
+        plt.title("ocl t=%.1f ms"%gpu_time)
+    if Ncpu > 0 and Ngpu > 0:
+        plt.subplot(133)
+        plot_data(data, 1e8*relerr)
+        plt.title("max rel err = %.3g"%max(abs(relerr)))
+        plt.colorbar()
+    plt.show()
+
+def rename(pars, **names):
+    newpars = pars.copy()
+    for new,old in names.items():
+        for variant in ("", "_pd", "_pd_n", "_pd_nsigma"):
+            if old+variant in newpars:
+                newpars[new+variant] = pars[old+variant]
+                del newpars[old+variant]
+    return newpars
+
+def rescale_sld(pars, names):
+    newpars = pars.copy()
+    for p in names:
+        newpars[p] *= 1e6
+    return newpars
+
+# ===========================================================================
+#
+MODELS = {}
+def model(name):
+    def gather_function(fn):
+        MODELS[name] = fn
+        return fn
+    return gather_function
+
+USAGE="""
+usage: compare model [Nopencl] [Nsasview]
+
+Compare the speed and value for a model between the SasView original and the
+OpenCL rewrite.
+
+* Nopencl is the number of times to run the OpenCL model (default=5)
+
+* Nsasview is the number of times to run the Sasview model (default=1)
+
+* model is the name of the model to compare:
+
+    %s
+"""
+
+def main():
+    if len(sys.argv) == 1:
+        models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' '))
+                               for k,v in sorted(MODELS.items()))
+        print(USAGE%models)
+        sys.exit(1)
+
+    cpuname, cpupars, gpuname, gpupars = MODELS[sys.argv[1]]()
+    Nopencl = int(sys.argv[2]) if len(sys.argv) > 2 else 5
+    Nsasview = int(sys.argv[3]) if len(sys.argv) > 3 else 1
+
+    compare(Nsasview, cpuname, cpupars, Nopencl, gpuname, gpupars)
+
+@model('cyl')
+def cylinder():
+    cpupars = dict(
+        scale=.003, background=.1,
+        sldCyl=.291e-6, sldSolv=5.77e-6,
+        radius=264.1, length=66.96,
+        cyl_theta=85, cyl_phi=0,
+        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
+        length_pd=0.1,length_pd_n=1, length_pd_nsigma=3,
+        cyl_theta_pd=45, cyl_theta_pd_n=50, cyl_theta_pd_nsigma=3,
+        cyl_phi_pd=0.1, cyl_phi_pd_n=5, cyl_phi_pd_nsigma=3,
+        )
+    cpuname = 'Cylinder'
+
+    gpupars = rename(cpupars, theta='cyl_theta', phi='cyl_phi', sld='sldCyl', solvent_sld='sldSolv')
+    gpupars = rescale_sld(gpupars, ['sld', 'solvent_sld'])
+    gpuname = 'cylinder'
+    return cpuname, cpupars, gpuname, gpupars
+
+
+@model('ell')
+def ellipse():
+    pars = dict(
+        scale=.027, background=4.9,
+        sldEll=.297e-6, sldSolv=5.773e-6,
+        radius_a=60, radius_b=180,
+        axis_theta=0, axis_phi=90,
+        radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3,
+        radius_b_pd=0.1, radius_b_pd_n=10, radius_b_pd_nsigma=3,
+        axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3,
+        axis_phi_pd=0.1, axis_phi_pd_n=6, axis_phi_pd_nsigma=3,
+        )
+
+    from Models.code_ellipse import GpuEllipse as gpumodel
+    model = sasview_model('Ellipsoid', **pars)
+
+    pars = rename(pars, theta='axis_theta', phi='axis_phi', sld='sldEll', solvent_sld='sldSolv')
+    pars = rescale_sld(pars, ['sld', 'solvent_sld'])
+    return model, gpumodel, pars
+
+
+@model('cscyl')
+def core_shell_cylinder(N=1):
+    pars = dict(
+        scale= 1.77881e-06, background=223.827,
+        core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6,
+        radius=325, thickness=25, length=34.2709,
+        axis_theta=90, axis_phi=0,
+        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
+        length_pd=0.1, length_pd_n=10, length_pd_nsigma=3,
+        thickness_pd=0.1, thickness_pd_n=5, thickness_pd_nsigma=3,
+        axis_theta_pd=15.8, axis_theta_pd_n=20, axis_theta_pd_nsigma=5,
+        axis_phi_pd=0.0008748, axis_phi_pd_n=5, axis_phi_pd_nsigma=3,
+        )
+
+    model = sasview_model('CoreShellCylinder', **pars)
+    from Models.code_coreshellcyl_f import GpuCoreShellCylinder as gpumodel
+
+    pars = rename(pars, theta='axis_theta', phi='axis_phi')
+    pars = rescale_sld(pars, ['core_sld', 'shell_sld', 'solvent_sld'])
+    return model, gpumodel, pars
+
+
+@model('ell3')
+def triaxial_ellipse(N=1):
+    pars = dict(
+        scale=0.08, background=5,
+        sldEll=7.105e-6, sldSolv=.291e-6,
+        axis_theta=0, axis_phi=0, axis_psi=0,
+        semi_axisA=15, semi_axisB=20, semi_axisC=500,
+        axis_theta_pd=20, axis_theta_pd_n=10, axis_theta_pd_nsigma=3,
+        axis_phi_pd=.1, axis_phi_pd_n=10, axis_phi_pd_nsigma=3,
+        axis_psi_pd=30, axis_psi_pd_n=5, axis_psi_pd_nsigma=3,
+        semi_axisA_pd=.1, semi_axisA_pd_n=5, semi_axisA_pd_nsigma=3,
+        semi_axisB_pd=.1, semi_axisB_pd_n=5, semi_axisB_pd_nsigma=3,
+        semi_axisC_pd=.1, semi_axisC_pd_n=5, semi_axisC_pd_nsigma=3,
+        )
+
+    model = sasview_model('TriaxialEllipsoid', **pars)
+    from Models.code_triaxialellipse import GpuTriEllipse as gpumodel
+
+    pars = rename(pars,
+                  theta='axis_theta', phi='axis_phi', psi='axis_psi',
+                  sld='sldEll', solvent_sld='sldSolv',
+                  radius_a='semi_axisA', radius_b='semi_axisB',
+                  radius_c='semi_axisC',
+                  )
+    pars = rescale_sld(pars, ['sld', 'solvent_sld'])
+    return model, gpumodel, pars
+
+@model('lam')
+def lamellar(N=1):
+    pars = dict(
+        scale=0.08, background=0.003,
+        sld_bi=5.38e-6,sld_sol=7.105e-6,
+        bi_thick=19.2946,
+        bi_thick_pd= 0.37765, bi_thick_pd_n=40, bi_thick_pd_nsigma=3,
+        )
+
+    model = sasview_model('Lamellar', **pars)
+    from Models.code_lamellar import GpuLamellar as gpumodel
+
+    pars = rename(pars, sld='sld_bi', solvent_sld='sld_sol', thickness='bi_thick')
+    pars = rescale_sld(pars, ['sld', 'solvent_sld'])
+    return model, gpumodel, pars
+
+@model('capcyl')
+def capped_cylinder(N=1):
+    pars = dict(
+        scale=.08, background=0,
+        sld_capcyl=1e-6, sld_solv=6.3e-6,
+        rad_cyl=20, rad_cap=40, len_cyl=400,
+        theta=0, phi=0,
+        rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_pd_nsigma=3,
+        rad_cap_pd=.1, rad_cap_pd_n=10, rad_cap_pd_nsigma=3,
+        len_cyl_pd=.1, len_cyl_pd_n=3, len_cyl_pd_nsigma=3,
+        theta_pd=.1, theta_pd_n=3, theta_pd_nsigma=3,
+        phi_pd=.1, phi_pd_n=3, phi_pd_nsigma=3,
+        )
+
+
+    model = sasview_model('CappedCylinder', **pars)
+    from Models.code_capcyl import GpuCapCylinder as gpumodel
+
+    pars = rename(pars,
+                  sld='sld_capcyl', solvent_sld='sld_solv',
+                  length='len_cyl', radius='rad_cyl',
+                  cap_radius='rad_cap')
+    pars = rescale_sld(pars, ['sld', 'solvent_sld'])
+    return model, gpumodel, pars
+
+
+if __name__ == "__main__":
+    main()
Index: sasmodels/alignment.py
===================================================================
--- sasmodels/alignment.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/alignment.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,54 @@
+"""
+GPU data alignment.
+
+Some web sites say that maximizing performance for OpenCL code requires
+aligning data on certain memory boundaries.
+
+:func:`data_alignment` queries all devices in the OpenCL context, returning
+the most restrictive alignment.
+
+:func:`align_data` aligns an existing array.
+
+:func:`align_empty` to create a new array of the correct alignment.
+
+Note:  This code is unused. So far, tests have not demonstrated any
+improvement from forcing correct alignment.  The tests should
+be repeated with arrays forced away from the target boundaries
+to decide whether it is really required.
+"""
+
+import numpy as np
+import pyopencl as cl
+
+def data_alignment(context):
+    """
+    Return the desired byte alignment across all devices.
+    """
+    # Note: rely on the fact that alignment will be 2^k
+    return max(d.min_data_type_align_size for d in context.devices)
+
+def align_empty(shape, dtype, alignment=128):
+    size = np.prod(shape)
+    dtype = np.dtype(dtype)
+    # allocate array with extra space for alignment
+    extra = alignment//dtype.itemsize - 1
+    result = np.empty(size+extra, dtype)
+    # build a view into allocated array which starts on a boundary
+    offset = (result.ctypes.data%alignment)//dtype.itemsize
+    view = np.reshape(result[offset:offset+size], shape)
+    return view
+
+def align_data(v, dtype, alignment=128):
+    # if v is contiguous, aligned, and of the correct type then just return v
+    view = align_empty(v.shape, dtype, alignment=alignment)
+    view[:] = v
+    return view
+
+def work_group_size(context, kernel):
+    """
+    Return the desired work group size for the context.
+    """
+    max(kernel.get_work_group_info(cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE, d)
+        for d in context.devices)
+
+
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/core.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,291 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys, os
+import datetime
+
+import numpy as np
+
+def load_model(modelname):
+    from os.path import abspath, join as joinpath, dirname
+    from sasmodels import gen
+    modelpath = abspath(joinpath(dirname(gen.__file__), 'models',
+                                 modelname+'.c'))
+    return gen.make(modelpath)
+
+def opencl_model(modelname, dtype="single"):
+    from sasmodels import gpu
+
+    source, meta, _ = 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)
+
+
+if sys.platform == 'darwin':
+    COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp"
+elif os.name == 'nt':
+    COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm"
+else:
+    COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm"
+DLL_PATH = "/tmp"
+def dll_path(meta):
+    from os.path import join as joinpath, split as splitpath, splitext
+    basename = splitext(splitpath(meta['filename'])[1])[0]
+    return joinpath(DLL_PATH, basename+'.so')
+
+def dll_model(modelname):
+    import os
+    from sasmodels import dll
+
+    source, meta, _ = load_model(modelname)
+    dllpath = dll_path(meta)
+    if not os.path.exists(dllpath) \
+            or (os.path.getmtime(dllpath) < os.path.getmtime(meta['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)
+
+TIC = None
+def tic():
+    global TIC
+    then = datetime.datetime.now()
+    TIC = lambda: (datetime.datetime.now()-then).total_seconds()
+    return TIC
+
+def toc():
+    return TIC()
+
+def load_data(filename):
+    from sans.dataloader.loader import Loader
+    loader = Loader()
+    data = loader.load(filename)
+    if data is None:
+        raise IOError("Data %r could not be loaded"%filename)
+    return data
+
+def fake_data2D(qx, qy=None):
+    from sans.dataloader.data_info import Data2D, Detector
+
+
+    if qy is None:
+        qy = qx
+    Qx,Qy = np.meshgrid(qx,qy)
+    Qx,Qy = Qx.flatten(), Qy.flatten()
+    Iq = 100*np.ones_like(Qx)
+    dIq = np.sqrt(Iq)
+    mask = np.ones(len(Iq), dtype='bool')
+
+    data = Data2D()
+    data.filename = "fake data"
+    data.qx_data = Qx
+    data.qy_data = Qy
+    data.data = Iq
+    data.err_data = dIq
+    data.mask = mask
+
+    # 5% dQ/Q resolution
+    data.dqx_data = 0.05*Qx
+    data.dqy_data = 0.05*Qy
+
+    detector = Detector()
+    detector.pixel_size.x = 5 # mm
+    detector.pixel_size.y = 5 # mm
+    detector.distance = 4 # m
+    data.detector.append(detector)
+    data.xbins = qx
+    data.ybins = qy
+    data.source.wavelength = 5 # angstroms
+    data.source.wavelength_unit = "A"
+    data.Q_unit = "1/A"
+    data.I_unit = "1/cm"
+    data.q_data = np.sqrt(Qx**2 + Qy**2)
+    data.xaxis("Q_x", "A^{-1}")
+    data.yaxis("Q_y", "A^{-1}")
+    data.zaxis("Intensity", r"\text{cm}^{-1}")
+    return data
+
+
+def set_beam_stop(data, radius, outer=None):
+    from sans.dataloader.manipulations import Ringcut
+    if hasattr(data, 'qx_data'):
+        data.mask = Ringcut(0, radius)(data)
+        if outer is not None:
+            data.mask += Ringcut(outer,np.inf)(data)
+    else:
+        data.mask = (data.x>=radius)
+        if outer is not None:
+            data.mask &= (data.x<outer)
+
+def set_half(data, half):
+    from sans.dataloader.manipulations import Boxcut
+    if half == 'right':
+        data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
+    if half == 'left':
+        data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
+
+def set_top(data, max):
+    from sans.dataloader.manipulations import Boxcut
+    data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data)
+
+def plot_data(data, iq, vmin=None, vmax=None):
+    from numpy.ma import masked_array
+    import matplotlib.pyplot as plt
+    img = masked_array(iq, data.mask)
+    xmin, xmax = min(data.qx_data), max(data.qx_data)
+    ymin, ymax = min(data.qy_data), max(data.qy_data)
+    plt.imshow(img.reshape(128,128),
+               interpolation='nearest', aspect=1, origin='upper',
+               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
+
+def plot_result2D(data, theory, view='linear'):
+    import matplotlib.pyplot as plt
+    from numpy.ma import masked_array, masked
+    #print "not a number",sum(np.isnan(data.data))
+    #data.data[data.data<0.05] = 0.5
+    mdata = masked_array(data.data, data.mask)
+    mdata[np.isnan(mdata)] = masked
+    if view is 'log':
+        mdata[mdata <= 0] = masked
+        mdata = np.log10(mdata)
+        mtheory = masked_array(np.log10(theory), mdata.mask)
+    else:
+        mtheory = masked_array(theory, mdata.mask)
+    mresid = masked_array((theory-data.data)/data.err_data, data.mask)
+    vmin = min(mdata.min(), mtheory.min())
+    vmax = max(mdata.max(), mtheory.max())
+    print np.exp(np.mean(mtheory)), np.std(mtheory),np.max(mtheory),np.min(mtheory)
+
+    plt.subplot(1, 3, 1)
+    plot_data(data, mdata, vmin=vmin, vmax=vmax)
+    plt.colorbar()
+    plt.subplot(1, 3, 2)
+    plot_data(data, mtheory, vmin=vmin, vmax=vmax)
+    plt.colorbar()
+    plt.subplot(1, 3, 3)
+    print abs(mresid).max()
+    plot_data(data, mresid)
+    plt.colorbar()
+
+
+def plot_result1D(data, theory, view='linear'):
+    import matplotlib.pyplot as plt
+    from numpy.ma import masked_array, masked
+    #print "not a number",sum(np.isnan(data.y))
+    #data.y[data.y<0.05] = 0.5
+    mdata = masked_array(data.y, data.mask)
+    mdata[np.isnan(mdata)] = masked
+    if view is 'log':
+        mdata[mdata <= 0] = masked
+    mtheory = masked_array(theory, mdata.mask)
+    mresid = masked_array((theory-data.y)/data.dy, mdata.mask)
+
+    plt.subplot(1,2,1)
+    plt.errorbar(data.x, mdata, yerr=data.dy)
+    plt.plot(data.x, mtheory, '-', hold=True)
+    plt.yscale(view)
+    plt.subplot(1, 2, 2)
+    plt.plot(data.x, mresid, 'x')
+    #plt.axhline(1, color='black', ls='--',lw=1, hold=True)
+    #plt.axhline(0, color='black', lw=1, hold=True)
+    #plt.axhline(-1, color='black', ls='--',lw=1, hold=True)
+
+
+class BumpsModel(object):
+    def __init__(self, data, model, cutoff=1e-5, **kw):
+        from bumps.names import Parameter
+        from . import gpu
+
+        # interpret data
+        self.is2D = hasattr(data,'qx_data')
+        self.data = data
+        if self.is2D:
+            self.index = (data.mask==0) & (~np.isnan(data.data))
+            self.iq = data.data[self.index]
+            self.diq = data.err_data[self.index]
+            self._theory = np.zeros_like(data.data)
+            q_vectors = [data.qx_data, data.qy_data]
+        else:
+            self.index = (data.mask==0) & (~np.isnan(data.y))
+            self.iq = data.y[self.index]
+            self.diq = data.dy[self.index]
+            self._theory = np.zeros_like(data.y)
+            q_vectors = [data.x]
+        #input = model.make_input(q_vectors)
+        input = model.make_input([v[self.index] for v in q_vectors])
+
+        # create model
+        self.fn = model(input, cutoff=cutoff)
+
+        # define bumps parameters
+        pars = []
+        for p in model.meta['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))
+                    pars.append(xname)
+        if kw:
+            raise TypeError("unexpected parameters: %s"%(", ".join(sorted(kw.keys()))))
+        self._parameter_names = pars
+        self.update()
+
+    def update(self):
+        self._cache = {}
+
+    def numpoints(self):
+        return len(self.iq)
+
+    def parameters(self):
+        return dict((k,getattr(self,k)) for k in self._parameter_names)
+
+    def theory(self):
+        if 'theory' not in self._cache:
+            pars = dict((k,getattr(self,k).value) for k in self._parameter_names)
+            #print pars
+            self._theory[self.index] = self.fn.eval(pars)
+            #self._theory[:] = self.fn.eval(pars)
+            self._cache['theory'] = self._theory
+        return self._cache['theory']
+
+    def residuals(self):
+        #if np.any(self.err ==0): print "zeros in err"
+        return (self.theory()[self.index]-self.iq)/self.diq
+
+    def nllf(self):
+        R = self.residuals()
+        #if np.any(np.isnan(R)): print "NaN in residuals"
+        return 0.5*np.sum(R**2)
+
+    def __call__(self):
+        return 2*self.nllf()/self.dof
+
+    def plot(self, view='log'):
+        if self.is2D:
+            plot_result2D(self.data, self.theory(), view=view)
+        else:
+            plot_result1D(self.data, self.theory(), view=view)
+
+    def save(self, basename):
+        pass
+
+def demo():
+    data = load_data('DEC07086.DAT')
+    set_beam_stop(data, 0.004)
+    plot_data(data)
+    import matplotlib.pyplot as plt; plt.show()
+
+if __name__ == "__main__":
+    demo()
Index: sasmodels/dll.py
===================================================================
--- sasmodels/dll.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/dll.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,124 @@
+"""
+C types wrapper for sasview models.
+"""
+
+import ctypes as ct
+from ctypes import c_void_p, c_int, c_double
+
+import numpy as np
+
+from . import gen
+
+from .gen import F32, F64
+
+IQ_ARGS = [c_void_p, c_void_p, c_int, c_void_p, c_double]
+IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int, c_void_p, c_double]
+
+class DllModel(object):
+    """
+    ctypes wrapper for a single model.
+
+    *source* and *meta* are the model source and interface as returned
+    from :func:`gen.make`.
+
+    *dtype* is the desired model precision.  Any numpy dtype for single
+    or double precision floats will do, such as 'f', 'float32' or 'single'
+    for single and 'd', 'float64' or 'double' for double.  Double precision
+    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)]
+
+
+        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] != ""]
+
+        # 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
+
+    def __call__(self, input, cutoff=1e-5):
+        kernel = self.Iqxy if input.is_2D else self.Iq
+        return DllKernel(kernel, self.meta, input, cutoff)
+
+    def make_input(self, q_vectors):
+        """
+        Make q input vectors available to the model.
+
+        This only needs to be done once for all models that operate on the
+        same input.  So for example, if you are adding two different models
+        together to compare to a data set, then only one model needs to
+        needs to call make_input, so long as the models have the same dtype.
+        """
+        return DllInput(q_vectors)
+
+
+class DllInput(object):
+    """
+    Make q data available to the gpu.
+
+    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
+    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
+    to get the best performance on OpenCL, which may involve shifting and
+    stretching the array to better match the memory architecture.  Additional
+    points will be evaluated with *q=1e-3*.
+
+    *dtype* is the data type for the q vectors. The data type should be
+    set to match that of the kernel, which is an attribute of
+    :class:`GpuProgram`.  Note that not all kernels support double
+    precision, so even if the program was created for double precision,
+    the *GpuProgram.dtype* may be single precision.
+
+    Call :meth:`release` when complete.  Even if not called directly, the
+    buffer will be released when the data object is freed.
+    """
+    def __init__(self, q_vectors):
+        self.nq = q_vectors[0].size
+        self.dtype = np.dtype('double')
+        self.is_2D = (len(q_vectors) == 2)
+        self.q_vectors = [np.ascontiguousarray(q,self.dtype) for q in q_vectors]
+        self.q_pointers = [q.ctypes.data for q in q_vectors]
+
+    def release(self):
+        self.q_vectors = []
+
+class DllKernel(object):
+    def __init__(self, kernel, meta, input, cutoff):
+        self.cutoff = cutoff
+        self.input = input
+        self.kernel = kernel
+        self.meta = meta
+
+        self.res = np.empty(input.nq, input.dtype)
+        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 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
+        #print pars
+        self.kernel(*pars)
+
+        return self.res
+
+    def release(self):
+        pass
+
+    def __del__(self):
+        self.release()
Index: sasmodels/gen.py
===================================================================
--- sasmodels/gen.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/gen.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,444 @@
+__all__ = ["make_opencl"]
+
+import os.path
+
+import numpy as np
+
+from .jsonutil import relaxed_loads
+
+F64 = np.dtype('float64')
+F32 = np.dtype('float32')
+
+# Conversion from units defined in the parameter table for each model
+# to units displayed in the sphinx documentation.
+RST_UNITS = {
+    "Ang": "|Ang|",
+    "1/Ang^2": "|Ang^-2|",
+    "1e-6/Ang^2": "|1e-6Ang^-2|",
+    "degrees": "degree",
+    "1/cm": "|cm^-1|",
+    "": "None",
+    }
+
+# Headers for the parameters tables in th sphinx documentation
+PARTABLE_HEADERS = [
+    "Parameter name",
+    "Units",
+    "Default value",
+    ]
+
+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.
+# This makes sure that the appropriate math constants are defined, and
+KERNEL_HEADER = """\
+// GENERATED CODE --- DO NOT EDIT ---
+// Code is produced by sasmodels.gen from sasmodels/models/MODEL.c
+
+#ifdef __OPENCL_VERSION__
+# define USE_OPENCL
+#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
+#  include <math.h>
+#  define REAL(x) (x)
+#  ifndef real
+#      define real double
+#  endif
+#  define global
+#  define local
+#  define constant const
+#  define kernel
+#  define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0)
+#else
+#  ifdef USE_SINCOS
+#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
+#  else
+#    define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0)
+#  endif
+#endif
+
+// Standard mathematical constants, prefixed with M_:
+//   E, LOG2E, LOG10E, LN2, LN10, PI, PI_2, PI_4, 1_PI, 2_PI,
+//   2_SQRTPI, SQRT2, SQRT1_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 REAL(3.141592653589793)
+#endif
+#ifndef M_PI_2
+#  define M_PI_2 REAL(1.570796326794897)
+#endif
+#ifndef M_PI_4
+#  define M_PI_4 REAL(0.7853981633974483)
+#endif
+
+// Non-standard pi/180, used for converting between degrees and radians
+#ifndef M_PI_180
+#  define M_PI_180 REAL(0.017453292519943295)
+#endif
+"""
+
+
+# The I(q) kernel and the I(qx, qy) kernel have one and two q parameters
+# respectively, so the template builder will need to do extra work to
+# declare, initialize and pass the q parameters.
+IQ_KERNEL = {
+    'fn': "Iq",
+    'q_par_decl': "global const real *q,",
+    'qinit': "const real qi = q[i];",
+    'qcall': "qi",
+    }
+
+IQXY_KERNEL = {
+    'fn': "Iqxy",
+    'q_par_decl': "global const real *qx,\n    global const real *qy,",
+    'qinit': "const real qxi = qx[i];\n    const real qyi = qy[i];",
+    'qcall': "qxi, qyi",
+    }
+
+# Generic kernel template for opencl/openmp.
+# This defines the opencl kernel that is available to the host.  The same
+# structure is used for Iq and Iqxy kernels, so extra flexibility is needed
+# for q parameters.  The polydispersity loop is built elsewhere and
+# substituted into this template.
+KERNEL_TEMPLATE = """\
+kernel void %(name)s(
+    %(q_par_decl)s
+    global real *result,
+#ifdef USE_OPENCL
+    global real *loops_g,
+#else
+    const int Nq,
+#endif
+    local real *loops,
+    const real cutoff,
+    %(par_decl)s
+    )
+{
+#ifdef USE_OPENCL
+  // copy loops info to local memory
+  event_t e = async_work_group_copy(loops, loops_g, (%(pd_length)s)*2, 0);
+  wait_group_events(1, &e);
+
+  int i = get_global_id(0);
+  int Nq = get_global_size(0);
+#endif
+
+#ifdef USE_OPENCL
+  if (i < Nq)
+#else
+  #pragma omp parallel for
+  for (int i=0; i < Nq; i++)
+#endif
+  {
+    %(qinit)s
+    real ret=REAL(0.0), norm=REAL(0.0);
+    real vol=REAL(0.0), norm_vol=REAL(0.0);
+%(loops)s
+    if (vol*norm_vol != REAL(0.0)) {
+      ret *= norm_vol/vol;
+    }
+    result[i] = scale*ret/norm+background;
+  }
+}
+"""
+
+# Polydispersity loop level.
+# This pulls the parameter value and weight from the looping vector in order
+# in preperation for a nested loop.
+LOOP_OPEN="""\
+for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) {
+  const real %(name)s = loops[2*(%(name)s_i%(offset)s)];
+  const real %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];"""
+
+# Polydispersity loop body.
+# This computes the weight, and if it is sufficient, calls the scattering
+# function and adds it to the total.  If there is a volume normalization,
+# it will also be added here.
+LOOP_BODY="""\
+const real weight = %(weight_product)s;
+if (weight > cutoff) {
+  ret += weight*%(fn)s(%(qcall)s, %(pcall)s);
+  norm += weight;
+  %(volume_norm)s
+}"""
+
+# 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
+# a normalized weight.
+VOLUME_NORM="""const real vol_weight = %(weight)s;
+  vol += vol_weight*volume(%(pars)s);
+  norm_vol += vol_weight;"""
+
+def indent(s, depth):
+    """
+    Indent a string of text with *depth* additional spaces on each line.
+    """
+    spaces = " "*depth
+    sep = "\n"+spaces
+    return spaces + sep.join(s.split("\n"))
+
+
+def make_kernel(meta, form):
+    """
+    Build a kernel call from metadata supplied by the user.
+
+    *meta* is the json object defined in the kernel file.
+
+    *form* is either "Iq" or "Iqxy".
+
+    This does not create a complete OpenCL kernel source, only the top
+    level kernel call with polydispersity and a call to the appropriate
+    Iq or Iqxy function.
+    """
+
+    # 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
+
+    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)
+        subst = { 'name': name, 'offset': offset }
+        loop_head.append(indent(LOOP_OPEN%subst, depth))
+        loop_end.insert(0, (" "*depth) + "}")
+        offset += '+N'+name
+        depth += 2
+
+    # The volume parameters in the inner loop are used to call the volume()
+    # function in the kernel, with the parameters defined in vol_pars and the
+    # weight product defined in weight.  If there are no volume parameters,
+    # then there will be no volume normalization.
+    if vol_pars:
+        subst = {
+            'weight': "*".join(p+"_w" for p in vol_pars),
+            'pars': ", ".join(vol_pars),
+            }
+        volume_norm = VOLUME_NORM%subst
+    else:
+        volume_norm = ""
+
+    # Define the inner loop function call
+    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),
+        }
+    loop_body = [indent(LOOP_BODY%subst, depth)]
+    loops = "\n".join(loop_head+loop_body+loop_end)
+
+    # declarations for non-pd followed by pd pars
+    # e.g.,
+    #     const real sld,
+    #     const int Nradius
+    fixed_par_decl = ",\n    ".join("const real %s"%p for p in fixed_pars)
+    pd_par_decl = ",\n    ".join("const int N%s"%p for p in pd_pars)
+    if fixed_par_decl and pd_par_decl:
+        par_decl = ",\n    ".join((fixed_par_decl, pd_par_decl))
+    elif fixed_par_decl:
+        par_decl = fixed_par_decl
+    else:
+        par_decl = pd_par_decl
+
+    # Finally, put the pieces together in the kernel.
+    subst = {
+        # kernel name is, e.g., cylinder_Iq
+        'name': "_".join((meta['name'], qpars['fn'])),
+        # to declare, e.g., global real q[],
+        'q_par_decl': qpars['q_par_decl'],
+        # to declare, e.g., real sld, int Nradius, int Nlength
+        'par_decl': par_decl,
+        # to copy global to local pd pars we need, e.g., Nradius+Nlength
+        'pd_length': "+".join('N'+p for p in pd_pars),
+        # the q initializers, e.g., real qi = q[i];
+        'qinit': qpars['qinit'],
+        # the actual polydispersity loop
+        'loops': loops,
+        }
+    kernel = KERNEL_TEMPLATE%subst
+    return kernel
+
+def make_partable(meta):
+    pars = meta['parameters']
+    column_widths = [
+        max(len(p[0]) for p in pars),
+        max(len(RST_UNITS[p[1]]) for p in pars),
+        PARTABLE_VALUE_WIDTH,
+        ]
+    column_widths = [max(w, len(h))
+                     for w,h in zip(column_widths, PARTABLE_HEADERS)]
+
+    sep = " ".join("="*w for w in column_widths)
+    lines = [
+        sep,
+        " ".join("%-*s"%(w,h) for w,h in zip(column_widths, PARTABLE_HEADERS)),
+        sep,
+        ]
+    for p in pars:
+        lines.append(" ".join([
+            "%-*s"%(column_widths[0],p[0]),
+            "%-*s"%(column_widths[1],RST_UNITS[p[1]]),
+            "%*g"%(column_widths[2],p[2]),
+            ]))
+    lines.append(sep)
+    return "\n".join(lines)
+
+def make_doc(kernelfile, meta, doc):
+    doc = doc%{'parameters': make_partable(meta)}
+    return doc
+
+def make_model(kernelfile, meta, source):
+    kernel_Iq = make_kernel(meta, "Iq")
+    kernel_Iqxy = make_kernel(meta, "Iqxy")
+    path = os.path.dirname(kernelfile)
+    extra = [open("%s/%s"%(path,f)).read() for f in meta['include']]
+    kernel = "\n\n".join([KERNEL_HEADER]+extra+[source, kernel_Iq, kernel_Iqxy])
+    return kernel
+
+def parse_file(kernelfile):
+    source = open(kernelfile).read()
+
+    # select parameters out of the source file
+    parts = source.split("PARAMETERS")
+    if len(parts) != 3:
+        raise ValueError("PARAMETERS block missing from %r"%kernelfile)
+    meta_source = parts[1].strip()
+    try:
+        meta = relaxed_loads(meta_source)
+    except:
+        print "in json text:"
+        print "\n".join("%2d: %s"%(i+1,s)
+                        for i,s in enumerate(meta_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())
+    elif len(parts) == 1:
+        raise ValueError("DOCUMENTATION block is missing from %r"%kernelfile)
+    else:
+        raise ValueError("DOCUMENTATION block incorrect from %r"%kernelfile)
+
+    return source, meta, doc
+
+def make(kernelfile):
+    """
+    Build an OpenCL function from the source in *kernelfile*.
+
+    The kernel file needs to define metadata about the parameters.  This
+    will be a JSON definition containing
+    """
+    #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
+
+
+def demo_time():
+    import datetime
+    tic = datetime.datetime.now()
+    toc = lambda: (datetime.datetime.now()-tic).total_seconds()
+    path = os.path.dirname("__file__")
+    doc, c = make_model(os.path.join(path, "models", "cylinder.c"))
+    print "time:",toc()
+
+def demo():
+    from os.path import join as joinpath, dirname
+    c, meta, doc = make_model(joinpath(dirname(__file__), "models", "cylinder.c"))
+    #print doc
+    #print c
+
+if __name__ == "__main__":
+    demo()
Index: sasmodels/gpu.py
===================================================================
--- sasmodels/gpu.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/gpu.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,265 @@
+"""
+GPU support through OpenCL
+
+There should be a single GPU environment running on the system.  This
+environment is constructed on the first call to :func:`env`, and the
+same environment is returned on each call.
+
+After retrieving the environment, the next step is to create the kernel.
+This is done with a call to :meth:`GpuEnvironment.make_kernel`, which
+returns the type of data used by the kernel.
+
+Next a :class:`GpuData` object should be created with the correct kind
+of data.  This data object can be used by multiple kernels, for example,
+if the target model is a weighted sum of multiple kernels.  The data
+should include any extra evaluation points required to compute the proper
+data smearing.  This need not match the square grid for 2D data if there
+is an index saying which q points are active.
+
+Together the GpuData, the program, and a device form a :class:`GpuKernel`.
+This kernel is used during fitting, receiving new sets of parameters and
+evaluating them.  The output value is stored in an output buffer on the
+devices, where it can be combined with other structure factors and form
+factors and have instrumental resolution effects applied.
+
+
+"""
+import warnings
+
+import numpy as np
+import pyopencl as cl
+from pyopencl import mem_flags as mf
+
+from . import gen
+
+from .gen import F32, F64
+
+F32_DEFS = """\
+#define REAL(x) (x##f)
+#define real float
+"""
+
+F64_DEFS = """\
+#pragma OPENCL EXTENSION cl_khr_fp64: enable
+#define REAL(x) (x)
+#define real double
+"""
+
+ENV = None
+def environment():
+    """
+    Returns a singleton :class:`GpuEnvironment`.
+
+    This provides an OpenCL context and one queue per device.
+    """
+    global ENV
+    if ENV is None:
+        ENV = GpuEnvironment()
+    return ENV
+
+def has_double(device):
+    """
+    Return true if device supports double precision.
+    """
+    return "cl_khr_fp64" in device.extensions
+
+
+def _stretch_input(vector, dtype, extra=1e-3, boundary=128):
+    """
+    Stretch an input vector to the correct boundary.
+
+    Performance on the kernels can drop by a factor of two or more if the
+    number of values to compute does not fall on a nice power of two
+    boundary.  A good choice for the boundary value is the
+    min_data_type_align_size property of the OpenCL device.  The usual
+    value of 128 gives a working size as a multiple of 32.  The trailing
+    additional vector elements are given a value of *extra*, and so
+    f(*extra*) will be computed for each of them.  The returned array
+    will thus be a subset of the computed array.
+    """
+    boundary // dtype.itemsize
+    remainder = vector.size%boundary
+    size = vector.size + (boundary - remainder if remainder != 0 else 0)
+    if size != vector.size:
+        vector = np.hstack((vector, [extra]*(size-vector.size)))
+    return np.ascontiguousarray(vector, dtype=dtype)
+
+
+def compile_model(context, source, dtype):
+    """
+    Build a model to run on the gpu.
+
+    Returns the compiled program and its type.  The returned type will
+    be float32 even if the desired type is float64 if any of the
+    devices in the context do not support the cl_khr_fp64 extension.
+    """
+    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
+
+    header = F64_DEFS if dtype == F64 else F32_DEFS
+    # Note: USE_SINCOS makes the intel cpu slower under opencl
+    if context.devices[0].type == cl.device_type.GPU:
+        header += "#define USE_SINCOS\n"
+    program  = cl.Program(context, header+source).build()
+    return program, dtype
+
+
+def make_result(self, size):
+    self.res = np.empty(size, dtype=self.dtype)
+    self.res_b = cl.Buffer(self.program.context, mf.READ_WRITE, self.res.nbytes)
+    return self.res, self.res_b
+
+
+# for now, this returns one device in the context
+# TODO: create a context that contains all devices on all platforms
+class GpuEnvironment(object):
+    """
+    GPU context, with possibly many devices, and one queue per device.
+    """
+    def __init__(self):
+        self.context = cl.create_some_context()
+        self.queues = [cl.CommandQueue(self.context, d)
+                       for d in self.context.devices]
+        self.boundary = max(d.min_data_type_align_size
+                            for d in self.context.devices)
+
+class GpuModel(object):
+    """
+    GPU wrapper for a single model.
+
+    *source* and *meta* are the model source and interface as returned
+    from :func:`gen.make`.
+
+    *dtype* is the desired model precision.  Any numpy dtype for single
+    or double precision floats will do, such as 'f', 'float32' or 'single'
+    for single and 'd', 'float64' or 'double' for double.  Double precision
+    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):
+        if self.dtype != input.dtype:
+            raise TypeError("data and kernel have different types")
+        kernel_name = gen.kernel_name(self.meta, input.is_2D)
+        kernel = getattr(self.program, kernel_name)
+        return GpuKernel(kernel, self.meta, input, cutoff)
+
+    def make_input(self, q_vectors):
+        """
+        Make q input vectors available to the model.
+
+        This only needs to be done once for all models that operate on the
+        same input.  So for example, if you are adding two different models
+        together to compare to a data set, then only one model needs to
+        needs to call make_input, so long as the models have the same dtype.
+        """
+        # Note: the weird interface, where make_input doesn't care which
+        # model calls it, allows us to ask the model to define the data
+        # and the caller does not need to know if it is opencl or ctypes.
+        # The returned data object is opaque.
+        return GpuInput(q_vectors, dtype=self.dtype)
+
+# TODO: check that we don't need a destructor for buffers which go out of scope
+class GpuInput(object):
+    """
+    Make q data available to the gpu.
+
+    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
+    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
+    to get the best performance on OpenCL, which may involve shifting and
+    stretching the array to better match the memory architecture.  Additional
+    points will be evaluated with *q=1e-3*.
+
+    *dtype* is the data type for the q vectors. The data type should be
+    set to match that of the kernel, which is an attribute of
+    :class:`GpuProgram`.  Note that not all kernels support double
+    precision, so even if the program was created for double precision,
+    the *GpuProgram.dtype* may be single precision.
+
+    Call :meth:`release` when complete.  Even if not called directly, the
+    buffer will be released when the data object is freed.
+    """
+    def __init__(self, q_vectors, dtype=F32):
+        env = environment()
+        self.nq = q_vectors[0].size
+        self.dtype = np.dtype(dtype)
+        self.is_2D = (len(q_vectors) == 2)
+        self.q_vectors = [
+            _stretch_input(q, self.dtype, boundary=env.boundary)
+            for q in q_vectors
+        ]
+        self.q_buffers = [
+            cl.Buffer(env.context,  mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
+            for q in self.q_vectors
+        ]
+        self.global_size = [self.q_vectors[0].size]
+
+    def release(self):
+        for b in self.q_buffers:
+            b.release()
+        self.q_buffers = []
+
+class GpuKernel(object):
+    def __init__(self, kernel, meta, input, cutoff):
+        env = environment()
+
+        self.cutoff = cutoff
+        self.input = input
+        self.kernel = kernel
+        self.meta = meta
+
+        # Inputs and outputs for each kernel call
+        self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE,
+                                  1024*input.dtype.itemsize)
+                        for _ in env.queues]
+        self.res_b = [cl.Buffer(env.context, mf.READ_WRITE,
+                                input.global_size[0]*input.dtype.itemsize)
+                      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):
+        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_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)
+        cl.enqueue_copy(queuei, self.res, res_bi)
+
+        return self.res
+
+    def release(self):
+        for b in self.loops_b:
+            b.release()
+        self.loops_b = []
+        for b in self.res_b:
+            b.release()
+        self.res_b = []
+
+    def __del__(self):
+        self.release()
Index: sasmodels/jsonutil.py
===================================================================
--- sasmodels/jsonutil.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/jsonutil.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,172 @@
+"""
+Read a relaxed JSON file.
+
+Relaxed JSON allows comments (introduced by // and going to the end
+of the line), optional quotes on key names in dictionaries and optional
+trailing commas at the ends of lists and dictionaries.  It also strips
+the leading characters up to the first '{'.  Multiline strings can be
+formatted using "\n\" at the end of each line.
+
+If the file contains e.g., "var _ = {...}", then it can be edited with
+a JavaScript aware editor such as VJET for eclipse, and it will be easier
+to locate errors and adjust formatting.
+"""
+import re
+import json
+from contextlib import contextmanager
+
+try:
+    from collections import OrderedDict
+except:
+    from ordered_dict import OrderedDict
+
+
+_LEADING_TEXT = re.compile(r'^.*?[{]', re.DOTALL)
+_LINE_CONTINUATION = re.compile(r'\\\s*\n')
+_TRAILING_COMMENT = re.compile(r'(?P<comment>\s*//.*?)\n')
+_MULTILINE_COMMENT = re.compile(r'(?P<comment>/\*.*?\*/)', re.DOTALL)
+_UNQUOTED_FIELDNAME = re.compile(r'(?P<prefix>[,{]\s*)(?P<key>[^\s,{}:"]+)(?P<tail>\s*:)')
+_TRAILING_COMMA = re.compile(r',(?P<tail>\s*[]}])')
+
+def relaxed_load(path, **kw):
+    return relaxed_loads(open(path).read(), **kw)
+
+def relaxed_loads(text, **kw):
+    """
+    Parse and return a relaxed JSON string.
+    """
+    ordered = kw.pop('ordered', False)
+    if ordered:  kw['object_pairs_hook'] = OrderedDict
+    # TODO: need a little state machine that performs the translation so that
+    # TODO: line and column numbers are preserved, and so that we can have
+    # TODO: http:// in a string (instead of it being treated like a comment). 
+    #print "== raw text\n", text
+    text = _LINE_CONTINUATION.sub('', text)
+    #print "== joined lines\n", text
+    text = _TRAILING_COMMENT.sub(r'\n', text)
+    text = _MULTILINE_COMMENT.sub(r'', text)
+    #print "== stripped comments\n", text
+    text = _LEADING_TEXT.sub('{', text)
+    #print "== trimmed text\n", text
+    text = _UNQUOTED_FIELDNAME.sub(r'\g<prefix>"\g<key>"\g<tail>', text)
+    #print "== quoted field names\n", text
+    text = _TRAILING_COMMA.sub(r'\g<tail>', text)
+    #print "== processed text\n", text
+    try:
+        obj = json.loads(text, object_hook=decode_dict_as_str, **kw)
+    except ValueError, e:
+        msg = [str(e)]
+        M = re.findall('line ([0-9]*) column ([0-9]*)', msg[0])
+        if M:
+            line,col = int(M[0][0]), int(M[0][1])
+            lines = text.split("\n")
+            if line>=2: msg.append(lines[line-2])
+            if line>=1: msg.append(lines[line-1])
+            msg.append(" "*(col-1) + "^")
+            if line<len(lines): msg.append(lines[line])
+        msg = "\n".join(msg)
+        raise e.__class__(msg)
+    return obj
+
+@contextmanager
+def float_format(formatstr='.15g'):
+    """
+    Allow the float format to be changed for a json encoding action.
+
+    This is a context manager, and should be used for example as::
+
+        >>> with float_format('.2g'):
+        >>>    print json.dumps(sqrt(2))
+        1.41 
+    """
+    formatter = json.encoder.FLOAT_REPR
+    json.encoder.FLOAT_REPR = lambda o: format(o, formatstr)
+    yield
+    json.encoder.FLOAT_REPR = formatter
+
+def numpy_encoder(o):
+    """
+    JSON encoder for numpy data.
+
+    To automatically convert numpy data to lists when writing a datastream
+    use json.dumps(object, default=numpy_json).
+    """
+    try:
+        return o.tolist()
+    except AttributeError:
+        raise TypeError
+
+
+def _decode_list(lst):
+    newlist = []
+    for i in lst:
+        if isinstance(i, unicode):
+            i = i.encode('utf-8')
+        elif isinstance(i, list):
+            i = _decode_list(i)
+        newlist.append(i)
+    return newlist
+
+def decode_dict_as_str(dct):
+    newdict = {}
+    for k, v in dct.iteritems():
+        if isinstance(k, unicode):
+            k = k.encode('utf-8')
+        if isinstance(v, unicode):
+            v = v.encode('utf-8')
+        elif isinstance(v, list):
+            v = _decode_list(v)
+        newdict[k] = v
+    return newdict
+
+
+def test():
+    """
+    Verify that the translation from pseudo-JSON to JSON works.
+    """
+    good = """\
+// This is a source definition with no errors
+var entry = {
+field : { // A comment about the field
+  "field" : "te\\
+x\\  
+t",
+  other$field : 56,
+  },
+/*
+multiline comment
+*/
+secondfield : {
+  content: ["string", "string"], /* a second comment */
+  content: [{name:"good", value:3, URL:"http:\\/\\/my.url.com"},]
+  },
+}
+"""
+    broken = """\
+// This is a source definition with a missing comma
+{
+field : { // A comment about the field
+  field : "te\\   
+x\\  
+t"
+  other$field : 56,
+  },
+/*
+multiline comment
+*/
+secondfield : {
+  content: ["string", "string"], /* a second comment */
+  content: [{name:"good", value:3},]
+  },
+}
+"""
+    result = relaxed_loads(good)
+    assert result['field']['field'] == "text"
+    assert result['field']['other$field'] == 56
+    assert result['secondfield']['content'][0]['name'] == 'good'
+    try: relaxed_loads(broken)
+    except ValueError, _: pass
+    else: raise Exception("No exception raised in broken")
+
+if __name__ == "__main__":
+    test()
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/models/cylinder.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,212 @@
+/* PARAMETERS
+{
+name: "cylinder",
+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" ],
+   [ "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" ],
+   [ "radius", "Ang",  20, [0, Infinity], "volume",
+     "Cylinder radius" ],
+   [ "length", "Ang",  400, [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
+.. _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 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)
+{
+    return M_PI*radius*radius*length;
+}
+
+
+real Iq(real q,
+    real sld,
+    real solvent_sld,
+    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 = (sld - solvent_sld) * volume(radius, length);
+    return REAL(1.0e-4) * form * s * s;
+}
+
+
+real Iqxy(real qx, real qy,
+    real sld,
+    real solvent_sld,
+    real radius,
+    real length,
+    real theta,
+    real phi)
+{
+    real sn, cn; // slots to hold sincos function output
+
+    // Compute angle alpha between q and the cylinder axis
+    SINCOS(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(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 = REAL(4.0)*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);
+    return REAL(1.0e-4) * form * s * s; // * correction;
+}
Index: sasmodels/models/lib/J1.c
===================================================================
--- sasmodels/models/lib/J1.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/models/lib/J1.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,38 @@
+real J1(real x);
+real J1(real x)
+{
+  const real ax = fabs(x);
+  if (ax < REAL(8.0)) {
+    const real y = x*x;
+    const real ans1 = x*(REAL(72362614232.0)
+              + y*(REAL(-7895059235.0)
+              + y*(REAL(242396853.1)
+              + y*(REAL(-2972611.439)
+              + y*(REAL(15704.48260)
+              + y*(REAL(-30.16036606)))))));
+    const real ans2 = REAL(144725228442.0)
+              + y*(REAL(2300535178.0)
+              + y*(REAL(18583304.74)
+              + y*(REAL(99447.43394)
+              + y*(REAL(376.9991397)
+              + y))));
+    return ans1/ans2;
+  } else {
+    const real y = REAL(64.0)/(ax*ax);
+    const real xx = ax - REAL(2.356194491);
+    const real ans1 = REAL(1.0)
+              + y*(REAL(0.183105e-2)
+              + y*(REAL(-0.3516396496e-4)
+              + y*(REAL(0.2457520174e-5)
+              + y*REAL(-0.240337019e-6))));
+    const real ans2 = REAL(0.04687499995)
+              + y*(REAL(-0.2002690873e-3)
+              + y*(REAL(0.8449199096e-5)
+              + y*(REAL(-0.88228987e-6)
+              + y*REAL(0.105787412e-6))));
+    real sn,cn;
+    SINCOS(xx, sn, cn);
+    const real ans = sqrt(REAL(0.636619772)/ax) * (cn*ans1 - (REAL(8.0)/ax)*sn*ans2);
+    return (x < REAL(0.0)) ? -ans : ans;
+  }
+}
Index: sasmodels/models/lib/cylkernel.c
===================================================================
--- sasmodels/models/lib/cylkernel.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/models/lib/cylkernel.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,16 @@
+// Compute the form factor for a cylinder
+// qq is the q-value for the calculation (1/A)
+// radius is the radius of the cylinder (A)
+// h is the HALF-LENGTH of the cylinder = L/2 (A)
+real CylKernel(real q, real radius, real h, real theta);
+real CylKernel(real q, real radius, real h, real theta)
+{
+    real sn, cn;
+    SINCOS(theta, sn, cn);
+    const real besarg = q*radius*sn;
+    const real siarg = q*h*cn;
+    // lim_{x->0} J1(x)/x = 1/2,  lim_{x->0} sin(x)/x = 1
+    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg);
+    const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg);
+    return REAL(4.0)*sin(theta)*bj*bj*si*si;
+}
Index: sasmodels/models/lib/gauss150.c
===================================================================
--- sasmodels/models/lib/gauss150.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/models/lib/gauss150.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,315 @@
+/*
+ *  GaussWeights.c
+ *  SANSAnalysis
+ *
+ *  Created by Andrew Jackson on 4/23/07.
+ *  Copyright 2007 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+// Gaussians
+constant real Gauss150Z[150]={
+  	REAL(-0.9998723404457334),
+  	REAL(-0.9993274305065947),
+  	REAL(-0.9983473449340834),
+  	REAL(-0.9969322929775997),
+  	REAL(-0.9950828645255290),
+  	REAL(-0.9927998590434373),
+  	REAL(-0.9900842691660192),
+  	REAL(-0.9869372772712794),
+  	REAL(-0.9833602541697529),
+  	REAL(-0.9793547582425894),
+  	REAL(-0.9749225346595943),
+  	REAL(-0.9700655145738374),
+  	REAL(-0.9647858142586956),
+  	REAL(-0.9590857341746905),
+  	REAL(-0.9529677579610971),
+  	REAL(-0.9464345513503147),
+  	REAL(-0.9394889610042837),
+  	REAL(-0.9321340132728527),
+  	REAL(-0.9243729128743136),
+  	REAL(-0.9162090414984952),
+  	REAL(-0.9076459563329236),
+  	REAL(-0.8986873885126239),
+  	REAL(-0.8893372414942055),
+  	REAL(-0.8795995893549102),
+  	REAL(-0.8694786750173527),
+  	REAL(-0.8589789084007133),
+  	REAL(-0.8481048644991847),
+  	REAL(-0.8368612813885015),
+  	REAL(-0.8252530581614230),
+  	REAL(-0.8132852527930605),
+  	REAL(-0.8009630799369827),
+  	REAL(-0.7882919086530552),
+  	REAL(-0.7752772600680049),
+  	REAL(-0.7619248049697269),
+  	REAL(-0.7482403613363824),
+  	REAL(-0.7342298918013638),
+  	REAL(-0.7198995010552305),
+  	REAL(-0.7052554331857488),
+  	REAL(-0.6903040689571928),
+  	REAL(-0.6750519230300931),
+  	REAL(-0.6595056411226444),
+  	REAL(-0.6436719971150083),
+  	REAL(-0.6275578900977726),
+  	REAL(-0.6111703413658551),
+  	REAL(-0.5945164913591590),
+  	REAL(-0.5776035965513142),
+  	REAL(-0.5604390262878617),
+  	REAL(-0.5430302595752546),
+  	REAL(-0.5253848818220803),
+  	REAL(-0.5075105815339176),
+  	REAL(-0.4894151469632753),
+  	REAL(-0.4711064627160663),
+  	REAL(-0.4525925063160997),
+  	REAL(-0.4338813447290861),
+  	REAL(-0.4149811308476706),
+  	REAL(-0.3959000999390257),
+  	REAL(-0.3766465660565522),
+  	REAL(-0.3572289184172501),
+  	REAL(-0.3376556177463400),
+  	REAL(-0.3179351925907259),
+  	REAL(-0.2980762356029071),
+  	REAL(-0.2780873997969574),
+  	REAL(-0.2579773947782034),
+  	REAL(-0.2377549829482451),
+  	REAL(-0.2174289756869712),
+  	REAL(-0.1970082295132342),
+  	REAL(-0.1765016422258567),
+  	REAL(-0.1559181490266516),
+  	REAL(-0.1352667186271445),
+  	REAL(-0.1145563493406956),
+  	REAL(-0.0937960651617229),
+  	REAL(-0.0729949118337358),
+  	REAL(-0.0521619529078925),
+  	REAL(-0.0313062657937972),
+  	REAL(-0.0104369378042598),
+  	REAL(0.0104369378042598),
+  	REAL(0.0313062657937972),
+  	REAL(0.0521619529078925),
+  	REAL(0.0729949118337358),
+  	REAL(0.0937960651617229),
+  	REAL(0.1145563493406956),
+  	REAL(0.1352667186271445),
+  	REAL(0.1559181490266516),
+  	REAL(0.1765016422258567),
+  	REAL(0.1970082295132342),
+  	REAL(0.2174289756869712),
+  	REAL(0.2377549829482451),
+  	REAL(0.2579773947782034),
+  	REAL(0.2780873997969574),
+  	REAL(0.2980762356029071),
+  	REAL(0.3179351925907259),
+  	REAL(0.3376556177463400),
+  	REAL(0.3572289184172501),
+  	REAL(0.3766465660565522),
+  	REAL(0.3959000999390257),
+  	REAL(0.4149811308476706),
+  	REAL(0.4338813447290861),
+  	REAL(0.4525925063160997),
+  	REAL(0.4711064627160663),
+  	REAL(0.4894151469632753),
+  	REAL(0.5075105815339176),
+  	REAL(0.5253848818220803),
+  	REAL(0.5430302595752546),
+  	REAL(0.5604390262878617),
+  	REAL(0.5776035965513142),
+  	REAL(0.5945164913591590),
+  	REAL(0.6111703413658551),
+  	REAL(0.6275578900977726),
+  	REAL(0.6436719971150083),
+  	REAL(0.6595056411226444),
+  	REAL(0.6750519230300931),
+  	REAL(0.6903040689571928),
+  	REAL(0.7052554331857488),
+  	REAL(0.7198995010552305),
+  	REAL(0.7342298918013638),
+  	REAL(0.7482403613363824),
+  	REAL(0.7619248049697269),
+  	REAL(0.7752772600680049),
+  	REAL(0.7882919086530552),
+  	REAL(0.8009630799369827),
+  	REAL(0.8132852527930605),
+  	REAL(0.8252530581614230),
+  	REAL(0.8368612813885015),
+  	REAL(0.8481048644991847),
+  	REAL(0.8589789084007133),
+  	REAL(0.8694786750173527),
+  	REAL(0.8795995893549102),
+  	REAL(0.8893372414942055),
+  	REAL(0.8986873885126239),
+  	REAL(0.9076459563329236),
+  	REAL(0.9162090414984952),
+  	REAL(0.9243729128743136),
+  	REAL(0.9321340132728527),
+  	REAL(0.9394889610042837),
+  	REAL(0.9464345513503147),
+  	REAL(0.9529677579610971),
+  	REAL(0.9590857341746905),
+  	REAL(0.9647858142586956),
+  	REAL(0.9700655145738374),
+  	REAL(0.9749225346595943),
+  	REAL(0.9793547582425894),
+  	REAL(0.9833602541697529),
+  	REAL(0.9869372772712794),
+  	REAL(0.9900842691660192),
+  	REAL(0.9927998590434373),
+  	REAL(0.9950828645255290),
+  	REAL(0.9969322929775997),
+  	REAL(0.9983473449340834),
+  	REAL(0.9993274305065947),
+  	REAL(0.9998723404457334)
+};
+
+constant real Gauss150Wt[150]={
+  	REAL(0.0003276086705538),
+  	REAL(0.0007624720924706),
+  	REAL(0.0011976474864367),
+  	REAL(0.0016323569986067),
+  	REAL(0.0020663664924131),
+  	REAL(0.0024994789888943),
+  	REAL(0.0029315036836558),
+  	REAL(0.0033622516236779),
+  	REAL(0.0037915348363451),
+  	REAL(0.0042191661429919),
+  	REAL(0.0046449591497966),
+  	REAL(0.0050687282939456),
+  	REAL(0.0054902889094487),
+  	REAL(0.0059094573005900),
+  	REAL(0.0063260508184704),
+  	REAL(0.0067398879387430),
+  	REAL(0.0071507883396855),
+  	REAL(0.0075585729801782),
+  	REAL(0.0079630641773633),
+  	REAL(0.0083640856838475),
+  	REAL(0.0087614627643580),
+  	REAL(0.0091550222717888),
+  	REAL(0.0095445927225849),
+  	REAL(0.0099300043714212),
+  	REAL(0.0103110892851360),
+  	REAL(0.0106876814158841),
+  	REAL(0.0110596166734735),
+  	REAL(0.0114267329968529),
+  	REAL(0.0117888704247183),
+  	REAL(0.0121458711652067),
+  	REAL(0.0124975796646449),
+  	REAL(0.0128438426753249),
+  	REAL(0.0131845093222756),
+  	REAL(0.0135194311690004),
+  	REAL(0.0138484622795371),
+  	REAL(0.0141714592928592),
+  	REAL(0.0144882814685445),
+  	REAL(0.0147987907597169),
+  	REAL(0.0151028518701744),
+  	REAL(0.0154003323133401),
+  	REAL(0.0156911024699895),
+  	REAL(0.0159750356447283),
+  	REAL(0.0162520081211971),
+  	REAL(0.0165218992159766),
+  	REAL(0.0167845913311726),
+  	REAL(0.0170399700056559),
+  	REAL(0.0172879239649355),
+  	REAL(0.0175283451696437),
+  	REAL(0.0177611288626114),
+  	REAL(0.0179861736145128),
+  	REAL(0.0182033813680609),
+  	REAL(0.0184126574807331),
+  	REAL(0.0186139107660094),
+  	REAL(0.0188070535331042),
+  	REAL(0.0189920016251754),
+  	REAL(0.0191686744559934),
+  	REAL(0.0193369950450545),
+  	REAL(0.0194968900511231),
+  	REAL(0.0196482898041878),
+  	REAL(0.0197911283358190),
+  	REAL(0.0199253434079123),
+  	REAL(0.0200508765398072),
+  	REAL(0.0201676730337687),
+  	REAL(0.0202756819988200),
+  	REAL(0.0203748563729175),
+  	REAL(0.0204651529434560),
+  	REAL(0.0205465323660984),
+  	REAL(0.0206189591819181),
+  	REAL(0.0206824018328499),
+  	REAL(0.0207368326754401),
+  	REAL(0.0207822279928917),
+  	REAL(0.0208185680053983),
+  	REAL(0.0208458368787627),
+  	REAL(0.0208640227312962),
+  	REAL(0.0208731176389954),
+  	REAL(0.0208731176389954),
+  	REAL(0.0208640227312962),
+  	REAL(0.0208458368787627),
+  	REAL(0.0208185680053983),
+  	REAL(0.0207822279928917),
+  	REAL(0.0207368326754401),
+  	REAL(0.0206824018328499),
+  	REAL(0.0206189591819181),
+  	REAL(0.0205465323660984),
+  	REAL(0.0204651529434560),
+  	REAL(0.0203748563729175),
+  	REAL(0.0202756819988200),
+  	REAL(0.0201676730337687),
+  	REAL(0.0200508765398072),
+  	REAL(0.0199253434079123),
+  	REAL(0.0197911283358190),
+  	REAL(0.0196482898041878),
+  	REAL(0.0194968900511231),
+  	REAL(0.0193369950450545),
+  	REAL(0.0191686744559934),
+  	REAL(0.0189920016251754),
+  	REAL(0.0188070535331042),
+  	REAL(0.0186139107660094),
+  	REAL(0.0184126574807331),
+  	REAL(0.0182033813680609),
+  	REAL(0.0179861736145128),
+  	REAL(0.0177611288626114),
+  	REAL(0.0175283451696437),
+  	REAL(0.0172879239649355),
+  	REAL(0.0170399700056559),
+  	REAL(0.0167845913311726),
+  	REAL(0.0165218992159766),
+  	REAL(0.0162520081211971),
+  	REAL(0.0159750356447283),
+  	REAL(0.0156911024699895),
+  	REAL(0.0154003323133401),
+  	REAL(0.0151028518701744),
+  	REAL(0.0147987907597169),
+  	REAL(0.0144882814685445),
+  	REAL(0.0141714592928592),
+  	REAL(0.0138484622795371),
+  	REAL(0.0135194311690004),
+  	REAL(0.0131845093222756),
+  	REAL(0.0128438426753249),
+  	REAL(0.0124975796646449),
+  	REAL(0.0121458711652067),
+  	REAL(0.0117888704247183),
+  	REAL(0.0114267329968529),
+  	REAL(0.0110596166734735),
+  	REAL(0.0106876814158841),
+  	REAL(0.0103110892851360),
+  	REAL(0.0099300043714212),
+  	REAL(0.0095445927225849),
+  	REAL(0.0091550222717888),
+  	REAL(0.0087614627643580),
+  	REAL(0.0083640856838475),
+  	REAL(0.0079630641773633),
+  	REAL(0.0075585729801782),
+  	REAL(0.0071507883396855),
+  	REAL(0.0067398879387430),
+  	REAL(0.0063260508184704),
+  	REAL(0.0059094573005900),
+  	REAL(0.0054902889094487),
+  	REAL(0.0050687282939456),
+  	REAL(0.0046449591497966),
+  	REAL(0.0042191661429919),
+  	REAL(0.0037915348363451),
+  	REAL(0.0033622516236779),
+  	REAL(0.0029315036836558),
+  	REAL(0.0024994789888943),
+  	REAL(0.0020663664924131),
+  	REAL(0.0016323569986067),
+  	REAL(0.0011976474864367),
+  	REAL(0.0007624720924706),
+  	REAL(0.0003276086705538)
+};
Index: sasmodels/models/lib/gauss20.c
===================================================================
--- sasmodels/models/lib/gauss20.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/models/lib/gauss20.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,55 @@
+/*
+ *  GaussWeights.c
+ *  SANSAnalysis
+ *
+ *  Created by Andrew Jackson on 4/23/07.
+ *  Copyright 2007 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+// Gaussians
+constant real Gauss20Wt[20]={
+	REAL(.0176140071391521),
+	REAL(.0406014298003869),
+	REAL(.0626720483341091),
+	REAL(.0832767415767047),
+	REAL(.10193011981724),
+	REAL(.118194531961518),
+	REAL(.131688638449177),
+	REAL(.142096109318382),
+	REAL(.149172986472604),
+	REAL(.152753387130726),
+	REAL(.152753387130726),
+	REAL(.149172986472604),
+	REAL(.142096109318382),
+	REAL(.131688638449177),
+	REAL(.118194531961518),
+	REAL(.10193011981724),
+	REAL(.0832767415767047),
+	REAL(.0626720483341091),
+	REAL(.0406014298003869),
+	REAL(.0176140071391521)
+};
+
+constant real Gauss20Z[20]={
+	REAL(-.993128599185095),
+	REAL(-.963971927277914),
+	REAL(-.912234428251326),
+	REAL(-.839116971822219),
+	REAL(-.746331906460151),
+	REAL(-.636053680726515),
+	REAL(-.510867001950827),
+	REAL(-.37370608871542),
+	REAL(-.227785851141645),
+	REAL(-.076526521133497),
+	REAL(.0765265211334973),
+	REAL(.227785851141645),
+	REAL(.37370608871542),
+	REAL(.510867001950827),
+	REAL(.636053680726515),
+	REAL(.746331906460151),
+	REAL(.839116971822219),
+	REAL(.912234428251326),
+	REAL(.963971927277914),
+	REAL(.993128599185095)
+};
Index: sasmodels/models/lib/gauss76.c
===================================================================
--- sasmodels/models/lib/gauss76.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/models/lib/gauss76.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,167 @@
+/*
+ *  GaussWeights.c
+ *  SANSAnalysis
+ *
+ *  Created by Andrew Jackson on 4/23/07.
+ *  Copyright 2007 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+// Gaussians
+constant real Gauss76Wt[76]={
+	REAL(.00126779163408536),		//0
+	REAL(.00294910295364247),
+	REAL(.00462793522803742),
+	REAL(.00629918049732845),
+	REAL(.00795984747723973),
+	REAL(.00960710541471375),
+	REAL(.0112381685696677),
+	REAL(.0128502838475101),
+	REAL(.0144407317482767),
+	REAL(.0160068299122486),
+	REAL(.0175459372914742),		//10
+	REAL(.0190554584671906),
+	REAL(.020532847967908),
+	REAL(.0219756145344162),
+	REAL(.0233813253070112),
+	REAL(.0247476099206597),
+	REAL(.026072164497986),
+	REAL(.0273527555318275),
+	REAL(.028587223650054),
+	REAL(.029773487255905),
+	REAL(.0309095460374916),		//20
+	REAL(.0319934843404216),
+	REAL(.0330234743977917),
+	REAL(.0339977794120564),
+	REAL(.0349147564835508),
+	REAL(.0357728593807139),
+	REAL(.0365706411473296),
+	REAL(.0373067565423816),
+	REAL(.0379799643084053),
+	REAL(.0385891292645067),
+	REAL(.0391332242205184),		//30
+	REAL(.0396113317090621),
+	REAL(.0400226455325968),
+	REAL(.040366472122844),
+	REAL(.0406422317102947),
+	REAL(.0408494593018285),
+	REAL(.040987805464794),
+	REAL(.0410570369162294),
+	REAL(.0410570369162294),
+	REAL(.040987805464794),
+	REAL(.0408494593018285),		//40
+	REAL(.0406422317102947),
+	REAL(.040366472122844),
+	REAL(.0400226455325968),
+	REAL(.0396113317090621),
+	REAL(.0391332242205184),
+	REAL(.0385891292645067),
+	REAL(.0379799643084053),
+	REAL(.0373067565423816),
+	REAL(.0365706411473296),
+	REAL(.0357728593807139),		//50
+	REAL(.0349147564835508),
+	REAL(.0339977794120564),
+	REAL(.0330234743977917),
+	REAL(.0319934843404216),
+	REAL(.0309095460374916),
+	REAL(.029773487255905),
+	REAL(.028587223650054),
+	REAL(.0273527555318275),
+	REAL(.026072164497986),
+	REAL(.0247476099206597),		//60
+	REAL(.0233813253070112),
+	REAL(.0219756145344162),
+	REAL(.020532847967908),
+	REAL(.0190554584671906),
+	REAL(.0175459372914742),
+	REAL(.0160068299122486),
+	REAL(.0144407317482767),
+	REAL(.0128502838475101),
+	REAL(.0112381685696677),
+	REAL(.00960710541471375),		//70
+	REAL(.00795984747723973),
+	REAL(.00629918049732845),
+	REAL(.00462793522803742),
+	REAL(.00294910295364247),
+	REAL(.00126779163408536)		//75 (indexed from 0)
+};
+
+constant real Gauss76Z[76]={
+	REAL(-.999505948362153),		//0
+	REAL(-.997397786355355),
+	REAL(-.993608772723527),
+	REAL(-.988144453359837),
+	REAL(-.981013938975656),
+	REAL(-.972229228520377),
+	REAL(-.961805126758768),
+	REAL(-.949759207710896),
+	REAL(-.936111781934811),
+	REAL(-.92088586125215),
+	REAL(-.904107119545567),		//10
+	REAL(-.885803849292083),
+	REAL(-.866006913771982),
+	REAL(-.844749694983342),
+	REAL(-.822068037328975),
+	REAL(-.7980001871612),
+	REAL(-.77258672828181),
+	REAL(-.74587051350361),
+	REAL(-.717896592387704),
+	REAL(-.688712135277641),
+	REAL(-.658366353758143),		//20
+	REAL(-.626910417672267),
+	REAL(-.594397368836793),
+	REAL(-.560882031601237),
+	REAL(-.526420920401243),
+	REAL(-.491072144462194),
+	REAL(-.454895309813726),
+	REAL(-.417951418780327),
+	REAL(-.380302767117504),
+	REAL(-.342012838966962),
+	REAL(-.303146199807908),		//30
+	REAL(-.263768387584994),
+	REAL(-.223945802196474),
+	REAL(-.183745593528914),
+	REAL(-.143235548227268),
+	REAL(-.102483975391227),
+	REAL(-.0615595913906112),
+	REAL(-.0205314039939986),
+	REAL(.0205314039939986),
+	REAL(.0615595913906112),
+	REAL(.102483975391227),			//40
+	REAL(.143235548227268),
+	REAL(.183745593528914),
+	REAL(.223945802196474),
+	REAL(.263768387584994),
+	REAL(.303146199807908),
+	REAL(.342012838966962),
+	REAL(.380302767117504),
+	REAL(.417951418780327),
+	REAL(.454895309813726),
+	REAL(.491072144462194),		//50
+	REAL(.526420920401243),
+	REAL(.560882031601237),
+	REAL(.594397368836793),
+	REAL(.626910417672267),
+	REAL(.658366353758143),
+	REAL(.688712135277641),
+	REAL(.717896592387704),
+	REAL(.74587051350361),
+	REAL(.77258672828181),
+	REAL(.7980001871612),	//60
+	REAL(.822068037328975),
+	REAL(.844749694983342),
+	REAL(.866006913771982),
+	REAL(.885803849292083),
+	REAL(.904107119545567),
+	REAL(.92088586125215),
+	REAL(.936111781934811),
+	REAL(.949759207710896),
+	REAL(.961805126758768),
+	REAL(.972229228520377),		//70
+	REAL(.981013938975656),
+	REAL(.988144453359837),
+	REAL(.993608772723527),
+	REAL(.997397786355355),
+	REAL(.999505948362153)		//75
+};
Index: sasmodels/weights.py
===================================================================
--- sasmodels/weights.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/weights.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,26 @@
+import numpy as np
+
+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
+
+    def get_pars(self):
+        return self.__dict__
+
+    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."""
+        npts, width, nsigmas = self.npts, self.width, self.nsigmas
+        sigma = width * center if relative else 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)]
+        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
+        return x, px
+
Index: setup_py2exe.py
===================================================================
--- setup_py2exe.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ setup_py2exe.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
@@ -0,0 +1,322 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2006-2010, University of Maryland
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/ or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+
+# Author: James Krycka
+
+"""
+This script uses py2exe to create dist\bumps.exe and dist\bumps_gui.exe for
+running the Bumps application in either CLI or GUI mode.
+
+These executables start the application and import the rest of the application
+code stored in library.zip.  The python interpreter and other required python
+packages and dlls are also placed in the zip file.  Additional resource files
+that are needed when Bumps is run are copied to the dist directory tree.  On
+completion, the contents of the dist directory tree can be used by the Inno
+Setup Compiler (via a separate script) to build a Windows installer/uninstaller
+for deployment of the Bumps application.  For testing purposes, bumps.exe or
+bumps_gui.exe can be run from the dist directory.
+"""
+
+import os
+import sys
+
+sys.dont_write_bytecode = True
+
+# Force build before continuing
+os.system('"%s" setup.py build' % sys.executable)
+
+# Remove the current directory from the python path
+here = os.path.abspath(os.path.dirname(__file__))
+sys.path = [p for p in sys.path if os.path.abspath(p) != here]
+
+import glob
+
+from distutils.core import setup
+from distutils.util import get_platform
+
+# Augment the setup interface with the py2exe command and make sure the py2exe
+# option is passed to setup.
+import py2exe
+
+if len(sys.argv) == 1:
+    sys.argv.append('py2exe')
+
+# Put the build lib on the start of the path.
+# For packages with binary extensions, need platform.  If it is a pure
+# script library, use an empty platform string.
+platform = '.%s-%s' % (get_platform(), sys.version[:3])
+#platform = ''
+build_lib = os.path.abspath('build/lib' + platform)
+sys.path.insert(0, build_lib)
+
+# print "\n".join(sys.path)
+
+import wx
+import matplotlib
+matplotlib.use('WXAgg')
+import periodictable
+
+# Retrieve the application version string.
+import bumps
+version = bumps.__version__
+from bumps.gui.resources import resources as gui_resources
+
+# A manifest is required to be included in a py2exe image (or accessible as a
+# file in the image directory) when wxPython is included so that the Windows XP
+# theme is used when rendering wx widgets.  The manifest must be matched to the
+# version of Python that is being used.
+#
+# Create a manifest for use with Python 2.5 on Windows XP or Vista.  It is
+# adapted from the Python manifest file (C:\Python25\pythonw.exe.manifest).
+
+manifest_for_python25 = """
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+<assembly xmlns="urn:schemas-microsoft-com:asm.v1" manifestVersion="1.0">
+<assemblyIdentity
+    version="1.0.0.0"
+    processorArchitecture="x86"
+    name="%(prog)s"
+    type="win32"
+/>
+<description>%(prog)s</description>
+<dependency>
+    <dependentAssembly>
+        <assemblyIdentity
+            type="win32"
+            name="Microsoft.Windows.Common-Controls"
+            version="6.0.0.0"
+            processorArchitecture="X86"
+            publicKeyToken="6595b64144ccf1df"
+            language="*"
+        />
+    </dependentAssembly>
+</dependency>
+</assembly>
+"""
+
+# Create a manifest for use with Python 2.6 or 2.7 on Windows XP or Vista.
+
+manifest_for_python26 = """
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+<assembly xmlns="urn:schemas-microsoft-com:asm.v1" manifestVersion="1.0">
+  <assemblyIdentity
+    version="5.0.0.0"
+    processorArchitecture="x86"
+    name="%(prog)s"
+    type="win32">
+  </assemblyIdentity>
+  <description>%(prog)s</description>
+  <trustInfo xmlns="urn:schemas-microsoft-com:asm.v3">
+    <security>
+      <requestedPrivileges>
+        <requestedExecutionLevel
+          level="asInvoker"
+          uiAccess="false">
+        </requestedExecutionLevel>
+      </requestedPrivileges>
+    </security>
+  </trustInfo>
+  <dependency>
+    <dependentAssembly>
+      <assemblyIdentity
+        type="win32"
+        name="Microsoft.VC90.CRT"
+        version="9.0.21022.8"
+        processorArchitecture="x86"
+        publicKeyToken="1fc8b3b9a1e18e3b">
+      </assemblyIdentity>
+    </dependentAssembly>
+  </dependency>
+  <dependency>
+    <dependentAssembly>
+      <assemblyIdentity
+        type="win32"
+        name="Microsoft.Windows.Common-Controls"
+        version="6.0.0.0"
+        processorArchitecture="x86"
+        publicKeyToken="6595b64144ccf1df"
+        language="*">
+      </assemblyIdentity>
+    </dependentAssembly>
+  </dependency>
+</assembly>
+"""
+
+# Select the appropriate manifest to use.
+if sys.version_info >= (3, 0) or sys.version_info < (2, 5):
+    print("*** This script only works with Python 2.5, 2.6, or 2.7.")
+    sys.exit()
+elif sys.version_info >= (2, 6):
+    manifest = manifest_for_python26
+elif sys.version_info >= (2, 5):
+    manifest = manifest_for_python25
+
+# Create a list of all files to include along side the executable being built
+# in the dist directory tree.  Each element of the data_files list is a tuple
+# consisting of a path (relative to dist\) and a list of files in that path.
+data_files = []
+
+# Add resource files that need to reside in the same directory as the image.
+data_files.append(('.', [os.path.join('.', 'LICENSE.txt')]))
+data_files.append(('.', [os.path.join('.', 'README.txt')]))
+data_files.append(('.', [os.path.join('.', 'bin', 'bumps_launch.bat')]))
+
+# Add application specific data files from the bumps\bumps-data folder.
+data_files += gui_resources.data_files()
+
+# Add data files from the matplotlib\mpl-data folder and its subfolders.
+# For matploblib prior to version 0.99 see the examples at the end of the file.
+data_files += matplotlib.get_py2exe_datafiles()
+
+# Add data files from the periodictable\xsf folder.
+data_files += periodictable.data_files()
+
+# Add example directories and their files.  An empty directory is ignored.
+# Note that Inno Setup will determine where these files will be placed such as
+# C:\My Documents\... instead of the installation folder.
+for path in glob.glob(os.path.join('examples', '*')):
+    if os.path.isdir(path):
+        for file in glob.glob(os.path.join(path, '*.*')):
+            data_files.append((path, [file]))
+    else:
+        data_files.append(('examples', [path]))
+
+for path in glob.glob(os.path.join('doc', 'examples', '*')):
+    if os.path.isdir(path):
+        for file in glob.glob(os.path.join(path, '*.*')):
+            data_files.append((path, [file]))
+    else:
+        data_files.append(('doc', [path]))
+
+# Add PDF documentation to the dist staging directory.
+pdf = os.path.join('doc', 'Bumps.pdf')
+if os.path.isfile(pdf):
+    data_files.append(('doc', [pdf]))
+else:
+    print("*** %s not found - building frozen image without it ***" % pdf)
+
+# Add the Microsoft Visual C++ 2008 redistributable kit if we are building with
+# Python 2.6 or 2.7.  This kit will be installed on the target system as part
+# of the installation process for the frozen image.  Note that the Python 2.5
+# interpreter requires msvcr71.dll which is included in the Python25 package,
+# however, Python 2.6 and 2.7 require the msvcr90.dll but they do not bundle it
+# with the Python26 or Python27 package.  Thus, for Python 2.6 and later, the
+# appropriate dll must be present on the target system at runtime.
+if sys.version_info >= (2, 6):
+    pypath = os.path.dirname(sys.executable)
+    data_files.append(('.', [os.path.join(pypath, 'vcredist_x86.exe')]))
+
+# Specify required packages to bundle in the executable image.
+packages = ['numpy', 'scipy', 'matplotlib', 'pytz', 'pyparsing',
+            'periodictable', 'bumps.names', 'dream'
+            ]
+
+# Specify files to include in the executable image.
+includes = []
+
+# Specify files to exclude from the executable image.
+# - We can safely exclude Tk/Tcl and Qt modules because our app uses wxPython.
+# - We do not use ssl services so they are omitted.
+# - We can safely exclude the TkAgg matplotlib backend because our app uses
+#   "matplotlib.use('WXAgg')" to override the default matplotlib configuration.
+# - On the web it is widely recommended to exclude certain lib*.dll modules
+#   but this does not seem necessary any more (but adding them does not hurt).
+# - Python25 requires mscvr71.dll, however, Win XP includes this file.
+# - Since we do not support Win 9x systems, w9xpopen.dll is not needed.
+# - For some reason cygwin1.dll gets included by default, but it is not needed.
+
+excludes = ['Tkinter', 'PyQt4', '_ssl', '_tkagg', 'numpy.distutils.test']
+
+dll_excludes = ['libgdk_pixbuf-2.0-0.dll',
+                'libgobject-2.0-0.dll',
+                'libgdk-win32-2.0-0.dll',
+                'tcl84.dll',
+                'tk84.dll',
+                'QtGui4.dll',
+                'QtCore4.dll',
+                'msvcr71.dll',
+                'msvcp90.dll',
+                'w9xpopen.exe',
+                'cygwin1.dll']
+
+
+class Target(object):
+
+    """This class stores metadata about the distribution in a dictionary."""
+
+    def __init__(self, **kw):
+        self.__dict__.update(kw)
+        self.version = version
+
+clientCLI = Target(
+    name='Bumps',
+    description='Bumps CLI application',
+    # module to run on application start
+    script=os.path.join('bin', 'bumps_cli.py'),
+    dest_base='bumps',  # file name part of the exe file to create
+    # also need to specify in data_files
+    icon_resources=[
+        (1, os.path.join('bumps', 'gui', 'resources', 'bumps.ico'))],
+    bitmap_resources=[],
+    other_resources=[(24, 1, manifest % dict(prog='Bumps'))])
+
+clientGUI = Target(
+    name='Bumps',
+    description='Bumps GUI application',
+    # module to run on application start
+    script=os.path.join('bin', 'bumps_gui.py'),
+    dest_base='bumps_gui',  # file name part of the exe file to create
+    # also need to specify in data_files
+    icon_resources=[
+        (1, os.path.join('bumps', 'gui', 'resources', 'bumps.ico'))],
+    bitmap_resources=[],
+    other_resources=[(24, 1, manifest % dict(prog='Bumps'))])
+
+# Now we do the work to create a standalone distribution using py2exe.
+#
+# When the application is run in console mode, a console window will be created
+# to receive any logging or error messages and the application will then create
+# a separate GUI application window.
+#
+# When the application is run in windows mode, it will create a GUI application
+# window and no console window will be provided.  Output to stderr will be
+# written to <app-image-name>.log.
+setup(
+    console=[clientCLI],
+    windows=[clientGUI],
+    options={'py2exe': {
+        'packages': packages,
+        'includes': includes,
+        'excludes': excludes,
+        'dll_excludes': dll_excludes,
+        'compressed': 1,   # standard compression
+        'optimize': 0,     # no byte-code optimization
+        'dist_dir': "dist",  # where to put py2exe results
+        'xref': False,     # display cross reference (as html doc)
+        'bundle_files': 1,  # bundle python25.dll in library.zip
+    }
+    },
+    # Since we are building two exe's, do not put the shared library in each
+    # of them.  Instead create a single, separate library.zip file.
+    # zipfile=None,               # bundle library.zip in exe
+    data_files=data_files           # list of files to copy to dist directory
+)
