Index: compare-new.py
===================================================================
--- compare-new.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ compare-new.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -7,5 +7,7 @@
 
 from sasmodels.bumps_model import BumpsModel, plot_data, tic
-from sasmodels.gen import opencl_model, dll_model
+from sasmodels import gpu, dll
+from sasmodels.convert import revert_model
+
 
 def sasview_model(modelname, **pars):
@@ -13,5 +15,5 @@
     Load a sasview model given the model name.
     """
-    modelname = modelname+"Model"
+    modelname = modelname
     sans = __import__('sans.models.'+modelname)
     ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None)
@@ -36,5 +38,5 @@
     sasmodels = __import__('sasmodels.models.'+modelname)
     module = getattr(sasmodels.models, modelname, None)
-    kernel = opencl_model(module, dtype=dtype)
+    kernel = gpu.load_model(module, dtype=dtype)
     return kernel
 
@@ -42,22 +44,67 @@
     sasmodels = __import__('sasmodels.models.'+modelname)
     module = getattr(sasmodels.models, modelname, None)
-    kernel = dll_model(module, dtype=dtype)
+    kernel = dll.load_model(module, dtype=dtype)
     return kernel
 
-
-def compare(Ncpu, cpuname, cpupars, Ngpu, gpuname, gpupars):
-
+def randomize(p, v):
+    """
+    Randomizing pars using sasview names.
+    Guessing at angles and slds.
+    """
+    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigmas','_pd_type')):
+        return v
+    if any(s in p for s in ('theta','phi','psi')):
+        return np.random.randint(0,90)
+    elif any(s in p for s in ('sld','Sld','SLD')):
+        return np.random.rand()*1e-5
+    else:
+        return np.random.rand()
+
+def parlist(pars):
+    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
+
+def suppress_pd(pars):
+    """
+    Suppress theta_pd for now until the normalization is resolved.
+
+    May also suppress complete polydispersity of the model to test
+    models more quickly.
+    """
+    for p in pars:
+        if p.endswith("_pd"): pars[p] = 0
+
+def compare(gpuname, gpupars, Ncpu, Ngpu, opts):
     #from sasmodels.core import load_data
     #data = load_data('December/DEC07098.DAT')
-    from sasmodels.core import empty_data1D
-    data = empty_data1D(np.logspace(-4, -1, 128))
-    #from sasmodels.core import empty_2D, set_beam_stop
-    #data = empty_data2D(np.linspace(-0.05, 0.05, 128))
-    #set_beam_stop(data, 0.004)
+    qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05)
+    if "-1d" in opts:
+        from sasmodels.bumps_model import empty_data1D
+        qmax = np.log10(qmax)
+        data = empty_data1D(np.logspace(qmax-3, qmax, 128))
+    else:
+        from sasmodels.bumps_model import empty_data2D, set_beam_stop
+        data = empty_data2D(np.linspace(-qmax, qmax, 128))
+        set_beam_stop(data, 0.004)
     is2D = hasattr(data, 'qx_data')
-
+    dtype = 'double' if '-double' in opts else 'single'
+    cutoff_opts = [s for s in opts if s.startswith('-cutoff')]
+    cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5
+
+
+    if '-random' in opts:
+        gpupars = dict((p,randomize(p,v)) for p,v in gpupars.items())
+    if '-pars' in opts: print "pars",parlist(gpupars)
+    if '-mono' in opts: suppress_pd(gpupars)
+
+    cpuname, cpupars = revert_model(gpuname, gpupars)
+
+    try:
+        gpumodel = load_opencl(gpuname, dtype=dtype)
+    except Exception,exc:
+        print exc
+        print "... trying again with single precision"
+        gpumodel = load_opencl(gpuname, dtype='single')
+    model = BumpsModel(data, gpumodel, cutoff=cutoff, **gpupars)
     if Ngpu > 0:
-        gpumodel = load_opencl(gpuname, dtype='single')
-        model = BumpsModel(data, gpumodel, **gpupars)
         toc = tic()
         for i in range(Ngpu):
@@ -66,10 +113,11 @@
             gpu = model.theory()
         gpu_time = toc()*1000./Ngpu
-        print "ocl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[~np.isnan(gpu)]))
+        print "ocl t=%.1f ms, intensity=%f"%(gpu_time, sum(gpu[~np.isnan(gpu)]))
         #print max(gpu), min(gpu)
 
-    if 0 and Ncpu > 0: # Hack to compare ctypes vs. opencl
-        dllmodel = load_dll(gpuname, dtype='double')
-        model = BumpsModel(data, dllmodel, **gpupars)
+    comp = None
+    if Ncpu > 0 and "-dll" in opts:
+        dllmodel = load_dll(gpuname, dtype=dtype)
+        model = BumpsModel(data, dllmodel, cutoff=cutoff, **gpupars)
         toc = tic()
         for i in range(Ncpu):
@@ -77,10 +125,10 @@
             cpu = model.theory()
         cpu_time = toc()*1000./Ncpu
-        print "dll t=%.1f ms"%cpu_time
+        comp = "dll"
 
     elif 0: # 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)
+        oldmodel = SasModel(data, oldgpu, dtype=dtype, **cpupars)
         toc = tic()
         for i in range(Ngpu):
@@ -88,5 +136,5 @@
             cpu = oldmodel.theory()
         cpu_time = toc()*1000./Ngpu
-        print "old t=%.1f ms"%cpu_time
+        comp = "old"
 
     elif Ncpu > 0:
@@ -100,14 +148,20 @@
                 cpu = cpumodel.evalDistribution(data.x)
         cpu_time = toc()*1000./Ncpu
-        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index]))
-
+        comp = 'sasview'
+        #print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index]))
+
+    if comp:
+        print "%s t=%.1f ms, intensity=%f"%(comp, cpu_time, sum(cpu[model.index]))
     if Ngpu > 0 and Ncpu > 0:
-        print "gpu/cpu", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu))
+        #print "speedup %.2g"%(cpu_time/gpu_time)
+        #print "max |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
+        resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu)
+        resid[model.index] = (gpu - cpu)[model.index]
+        relerr[model.index] = resid[model.index]/cpu[model.index]
+        print "max(|ocl-%s|)"%comp, max(abs(resid[model.index]))
+        print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[model.index]))
+
+    if '-noplot' in opts: return
 
     import matplotlib.pyplot as plt
@@ -115,5 +169,5 @@
         if Ngpu > 0: plt.subplot(131)
         plot_data(data, cpu, scale='log')
-        plt.title("omp t=%.1f ms"%cpu_time)
+        plt.title("%s t=%.1f ms"%(comp,cpu_time))
     if Ngpu > 0:
         if Ncpu > 0: plt.subplot(132)
@@ -122,26 +176,81 @@
     if Ncpu > 0 and Ngpu > 0:
         plt.subplot(133)
-        plot_data(data, 1e8*relerr, scale='linear')
-        plt.title("max rel err = %.3g"%max(abs(relerr)))
+        err = resid if '-abs' in opts else relerr
+        plot_data(data, err, scale='linear')
+        plt.title("max rel err = %.3g"%max(abs(err[model.index])))
         if is2D: 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
-
 # ===========================================================================
 #
+USAGE="""
+usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val]
+
+Compare the speed and value for a model between the SasView original and the
+OpenCL rewrite.
+
+model is the name of the model to compare (see below).
+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)
+
+Options (* for default):
+
+    -plot*/-noplot plots or suppress the plot of the model
+    -single/-double uses double precision for comparison
+    -lowq/-midq/-highq use q values up to 0.05, 0.2 or 1.0
+    -2d/-1d uses 1d or 2d random data
+    -preset/-random randomizes the parameters
+    -poly/-mono force monodisperse/polydisperse
+    -sasview/-dll whether cpu is tested using sasview or dll
+    -cutoff=1e-5/value cutoff for including a point in polydispersity
+    -nopars/-pars prints the parameter set or not
+    -rel/-abs plot relative or absolute error
+
+Key=value pairs allow you to set specific values to any of the model
+parameters.
+
+Available models:
+
+    %s
+"""
+
+def main():
+    opts = [arg for arg in sys.argv[1:] if arg.startswith('-')]
+    args = [arg for arg in sys.argv[1:] if not arg.startswith('-')]
+    models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' '))
+                           for k,v in sorted(MODELS.items()))
+    if len(args) == 0:
+        print(USAGE%models)
+        sys.exit(1)
+    if args[0] not in MODELS:
+        print "Model %r not available. Use one of:\n    %s"%(args[0],models)
+        sys.exit(1)
+
+    name, pars = MODELS[args[0]]()
+    Nopencl = int(args[1]) if len(args) > 1 else 5
+    Nsasview = int(args[2]) if len(args) > 3 else 1
+
+    # Fill in default polydispersity parameters
+    pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd'))
+    for p in pds:
+        if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3
+        if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian"
+
+    # Fill in parameters given on the command line
+    for arg in args[3:]:
+        k,v = arg.split('=')
+        if k not in pars:
+            # extract base name without distribution
+            # style may be either a.d or a_pd_d
+            s = set((p.split('_pd')[0]).split('.')[0] for p in pars)
+            print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s)))
+            sys.exit(1)
+        pars[k] = float(v) if not v.endswith('type') else v
+
+    compare(name, pars, Nsasview, Nopencl, opts)
+
+# ===========================================================================
+#
+
 MODELS = {}
 def model(name):
@@ -151,163 +260,102 @@
     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
+    pars = dict(
+        scale=1, background=0,
+        sld=6e-6, solvent_sld=1e-6,
+        #radius=5, length=20,
+        radius=260, length=290,
+        theta=30, phi=15,
+        radius_pd=.2, radius_pd_n=1,
+        length_pd=.2,length_pd_n=1,
+        theta_pd=15, theta_pd_n=25,
+        phi_pd=15, phi_pd_n=1,
+        )
+    return 'cylinder', pars
+
+@model('capcyl')
+def capped_cylinder():
+    pars = dict(
+        scale=1, background=0,
+        sld=6e-6, solvent_sld=1e-6,
+        radius=260, cap_radius=80000, length=290,
+        theta=30, phi=15,
+        radius_pd=.2, radius_pd_n=1,
+        cap_radius_pd=.2, cap_radius_pd_n=1,
+        length_pd=.2, length_pd_n=1,
+        theta_pd=15, theta_pd_n=25,
+        phi_pd=15, phi_pd_n=1,
+        )
+    return 'capped_cylinder', pars
+
+
+@model('cscyl')
+def core_shell_cylinder():
+    pars = dict(
+        scale=1, background=0,
+        core_sld=6e-6, shell_sld=8e-6, solvent_sld=1e-6,
+        radius=325, thickness=25, length=34.2709,
+        theta=90, phi=0,
+        radius_pd=0.1, radius_pd_n=10,
+        length_pd=0.1, length_pd_n=5,
+        thickness_pd=0.1, thickness_pd_n=5,
+        theta_pd=15.8, theta_pd_n=5,
+        phi_pd=0.0008748, phi_pd_n=5,
+        )
+    return 'core_shell_cylinder', pars
 
 
 @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
+def ellipsoid():
+    pars = dict(
+        scale=1, background=0,
+        sld=6e-6, solvent_sld=1e-6,
+        rpolar=50, requatorial=30,
+        theta=0, phi=0,
+        rpolar_pd=0.3, rpolar_pd_n=10,
+        requatorial_pd=0, requatorial_pd_n=10,
+        theta_pd=0, theta_pd_n=45,
+        phi_pd=0, phi_pd_n=45,
+        )
+    return 'ellipsoid', 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
+def triaxial_ellipsoid():
+    pars = dict(
+        scale=1, background=0,
+        sld=6e-6, solvent_sld=1e-6,
+        theta=0, phi=0, psi=0,
+        req_minor=25, req_major=36, rpolar=50,
+        theta_pd=0, theta_pd_n=5,
+        phi_pd=0, phi_pd_n=5,
+        psi_pd=0, psi_pd_n=5,
+        req_minor_pd=0, req_minor_pd_n=5,
+        req_major_pd=0, req_major_pd_n=5,
+        rpolar_pd=.3, rpolar_pd_n=25,
+        )
+    return 'triaxial_ellipsoid', pars
+
+@model('sph')
+def sphere():
+    pars = dict(
+        scale=1, background=0,
+        sld=6e-6, solvent_sld=1e-6,
+        radius=120,
+        radius_pd=.3, radius_pd_n=5,
+        )
+    return 'sphere', 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
-
+def lamellar():
+    pars = dict(
+        scale=1, background=0,
+        sld=6e-6, solvent_sld=1e-6,
+        thickness=40,
+        thickness_pd= 0.3, thickness_pd_n=40,
+        )
+    return 'lamellar', pars
 
 if __name__ == "__main__":
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/bumps_model.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -140,11 +140,15 @@
     import matplotlib.pyplot as plt
     if hasattr(data, 'qx_data'):
-        img = masked_array(iq, data.mask)
+        iq = iq[:]
+        valid = np.isfinite(iq)
         if scale == 'log':
-            img[(img <= 0) | ~np.isfinite(img)] = masked
-            img = np.log10(img)
+            valid[valid] = (iq[valid] > 0)
+            iq[valid] = np.log10(iq[valid])
+        iq[~valid|data.mask] = 0
+        #plottable = iq
+        plottable = masked_array(iq, ~valid|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),
+        plt.imshow(plottable.reshape(128,128),
                    interpolation='nearest', aspect=1, origin='upper',
                    extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
@@ -154,5 +158,6 @@
             plt.plot(data.x[idx], iq[idx])
         else:
-            idx = np.isfinite(iq) & (iq>0)
+            idx = np.isfinite(iq)
+            idx[idx] = (iq[idx]>0)
             plt.loglog(data.x[idx], iq[idx])
 
@@ -228,4 +233,5 @@
     def __init__(self, data, model, cutoff=1e-5, **kw):
         from bumps.names import Parameter
+        partype = model.info['partype']
 
         # interpret data
@@ -236,5 +242,8 @@
             self.diq = data.err_data[self.index]
             self._theory = np.zeros_like(data.data)
-            q_vectors = [data.qx_data, data.qy_data]
+            if not partype['orientation'] and not partype['magnetic']:
+                q_vectors = [np.sqrt(data.qx_data**2+data.qy_data**2)]
+            else:
+                q_vectors = [data.qx_data, data.qy_data]
         else:
             self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y)
@@ -258,5 +267,5 @@
             setattr(self, name, Parameter.default(value, name=name, limits=limits))
             pars.append(name)
-        for name in model.info['partype']['pd-2d']:
+        for name in partype['pd-2d']:
             for xpart,xdefault,xlimits in [
                     ('_pd', 0, limits),
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/convert.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,108 @@
+"""
+Convert models to and from sasview.
+"""
+PARAMETER_MAP = {
+    'CylinderModel': dict(
+        name='cylinder',
+        cyl_theta='theta', cyl_phi='phi',
+        sldCyl='sld', sldSolv='solvent_sld',
+        ),
+    'EllipsoidModel': dict(
+        name='ellipsoid',
+        axis_theta='theta', axis_phi='phi',
+        sldEll='sld', sldSolv='solvent_sld',
+        radius_a='rpolar', radius_b='requatorial',
+        ),
+    'CoreShellCylinderModel': dict(
+        name='core_shell_cylinder',
+        axis_theta='theta', axis_phi='phi',
+        ),
+    'TriaxialEllipsoidModel': dict(
+        name='triaxial_ellipsoid',
+        axis_theta='theta', axis_phi='phi', axis_psi='psi',
+        sldEll='sld', sldSolv='solvent_sld',
+        semi_axisA='req_minor', semi_axisB='req_major', semi_axisC='rpolar',
+        ),
+    'LamellarModel': dict(
+        name='lamellar',
+        sld_bi='sld', sld_sol='solvent_sld',
+        bi_thick='thickness',
+        ),
+    'CappedCylinderModel': dict(
+        name='capped_cylinder',
+        sld_capcyl='sld', sld_solv='solvent_sld',
+        len_cyl='length', rad_cyl='radius', rad_cap='cap_radius',
+        ),
+    'SphereModel': dict(
+        name='sphere',
+        sldSph='sld', sldSolv='solvent_sld',
+        #radius='radius',  # DO NOT LIST IDENTICAL PARAMETERS!
+        ),
+    }
+
+def _reverse_map():
+    retval = {}
+    for old_name,old_map in PARAMETER_MAP.items():
+        new_name = old_map['name']
+        new_map = dict((v,k) for k,v in old_map.items() if k != 'name')
+        new_map['name'] = old_name
+        retval[new_name] = new_map
+    return retval
+REVERSE_MAP = _reverse_map()
+del _reverse_map
+
+
+def _rename_pars(pars, mapping):
+    """
+    Rename the parameters and any associated polydispersity attributes.
+    """
+    newpars = pars.copy()
+    for old,new in mapping.items():
+        # Bumps style parameter names
+        for variant in ("", "_pd", "_pd_n", "_pd_nsigma", "_pd_type"):
+            if old+variant in newpars:
+                newpars[new+variant] = pars[old+variant]
+                del newpars[old+variant]
+        # Sasview style parameter names
+        for variant in (".width", ".nsigmas", ".type", ".npts"):
+            if old+variant in newpars:
+                newpars[new+variant] = pars[old+variant]
+                del newpars[old+variant]
+    return newpars
+
+def _rescale_sld(pars):
+    """
+    rescale all sld parameters in the new model definition by 1e6 so the
+    numbers are nicer.  Relies on the fact that all sld parameters in the
+    new model definition end with sld.
+    """
+    return dict((p, (v*1e6 if p.endswith('sld') else v))
+                for p,v in pars.items())
+
+def convert_model(name, pars):
+    """
+    Convert model from old style parameter names to new style.
+    """
+    mapping = PARAMETER_MAP[name]
+    newname = mapping['name']
+    newpars = _rescale_sld(_rename_pars(pars, mapping))
+    return newname, newpars
+
+def _unscale_sld(pars):
+    """
+    rescale all sld parameters in the new model definition by 1e6 so the
+    numbers are nicer.  Relies on the fact that all sld parameters in the
+    new model definition end with sld.
+    """
+    return dict((p, (v*1e-6 if p.endswith('sld') else v))
+                for p,v in pars.items())
+
+def revert_model(name, pars):
+    """
+    Convert model from new style parameter names to old style.
+    """
+    mapping = REVERSE_MAP[name]
+    oldname = mapping['name']
+    oldpars = _rename_pars(_unscale_sld(pars), mapping)
+    return oldname, oldpars
+
Index: sasmodels/dll.py
===================================================================
--- sasmodels/dll.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/dll.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -2,5 +2,6 @@
 C types wrapper for sasview models.
 """
-
+import sys
+import os
 import ctypes as ct
 from ctypes import c_void_p, c_int, c_double
@@ -11,4 +12,53 @@
 
 from .gen import F32, F64
+# Compiler platform details
+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(info):
+    """
+    Path to the compiled model defined by *info*.
+    """
+    from os.path import join as joinpath, split as splitpath, splitext
+    basename = splitext(splitpath(info['filename'])[1])[0]
+    return joinpath(DLL_PATH, basename+'.so')
+
+
+def load_model(kernel_module, dtype=None):
+    """
+    Load the compiled model defined by *kernel_module*.
+
+    Recompile if any files are newer than the model file.
+
+    *dtype* is ignored.  Compiled files are always double.
+
+    The DLL is not loaded until the kernel is called so models an
+    be defined without using too many resources.
+    """
+    import tempfile
+
+    source, info = gen.make(kernel_module)
+    source_files = gen.sources(info) + [info['filename']]
+    newest = max(os.path.getmtime(f) for f in source_files)
+    dllpath = dll_path(info)
+    if not os.path.exists(dllpath) or os.path.getmtime(dllpath)<newest:
+        # Replace with a proper temp file
+        fid, filename = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name'])
+        os.fdopen(fid,"w").write(source)
+        status = os.system(COMPILE%(filename, dllpath))
+        if status != 0:
+            print "compile failed.  File is in %r"%filename
+        else:
+            ## uncomment the following to keep the generated c file
+            #os.unlink(filename); print "saving compiled file in %r"%filename
+            pass
+    return DllModel(dllpath, info)
+
 
 IQ_ARGS = [c_void_p, c_void_p, c_int, c_void_p, c_double]
@@ -124,4 +174,5 @@
     """
     def __init__(self, kernel, info, input):
+        self.info = info
         self.input = input
         self.kernel = kernel
Index: sasmodels/gen.py
===================================================================
--- sasmodels/gen.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/gen.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -28,26 +28,40 @@
 
 *Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99
-functions written for OpenCL.  Floating point values should be
-declared as *real*.  Depending on how the function is called, a macro
-will replace *real* with *float* or *double*.  Unfortunately, MacOSX
-is picky about floating point constants, which should be defined with
-value + 'f' if they are of type *float* or just a bare value if they
-are of type *double*.  The solution is a macro *REAL(value)* which
-adds the 'f' if compiling for single precision floating point.  This
-does make the code ugly, and may someday be replaced by a clever
-regular expression which does the same job.  OpenCL has a *sincos*
-function which can improve performance when both the *sin* and *cos*
-values are needed for a particular argument.  Since this function
-does not exist in C-99, all use of *sincos* should be replaced by the
-macro *SINCOS(value,sn,cn)* where *sn* and *cn* are previously declared
-*real* values.  *value* may be an expression.  When compiled for systems
-without OpenCL, *SINCOS* will be replaced by *sin* and *cos* calls.  All
-functions need prototype declarations even if the are defined before they
-are used -- another present from MacOSX.  OpenCL does not support
-*#include* preprocessor directives; instead the includes must be listed
-in the kernel metadata, with functions defined before they are used.
-The included files should be listed using relative path to the kernel
-source file, or if using one of the standard models, relative to the
-sasmodels source files.
+functions written for OpenCL.  All functions need prototype declarations
+even if the are defined before they are used.  OpenCL does not support
+*#include* preprocessor directives, so instead the list of includes needs
+to be given as part of the metadata in the kernel module definition.
+The included files should be listed using a path relative to the kernel
+module, or if using "lib/file.c" if it is one of the standard includes
+provided with the sasmodels source.  The includes need to be listed in
+order so that functions are defined before they are used.
+
+Floating point values should be declared as *real*.  Depending on how the
+function is called, a macro will replace *real* with *float* or *double*.
+Unfortunately, MacOSX is picky about floating point constants, which
+should be defined with value + 'f' if they are of type *float* or just
+a bare value if they are of type *double*.  The solution is a macro
+*REAL(value)* which adds the 'f' if compiling for single precision
+floating point.  [Note: we could write a clever regular expression
+which automatically detects real-valued constants.  If we wanted to
+make our code more C-like, we could define variables with double but
+replace them with float before compiling for single precision.]
+
+OpenCL has a *sincos* function which can improve performance when both
+the *sin* and *cos* values are needed for a particular argument.  Since
+this function does not exist in C99, all use of *sincos* should be
+replaced by the macro *SINCOS(value,sn,cn)* where *sn* and *cn* are
+previously declared *real* values.  *value* may be an expression.  When
+compiled for systems without OpenCL, *SINCOS* will be replaced by
+*sin* and *cos* calls.
+
+If the input parameters are invalid, the scattering calculator should
+return a negative number. Particularly with polydispersity, there are
+some sets of shape parameters which lead to nonsensical forms, such
+as a capped cylinder where the cap radius is smaller than the
+cylinder radius.  The polydispersity calculation will ignore these points,
+effectively chopping the parameter weight distributions at the boundary
+of the infeasible region.  The resulting scattering will be set to
+background.  This will work correctly even when polydispersity is off.
 
 *ER* and *VR* are python functions which operate on parameter vectors.
@@ -112,4 +126,8 @@
     *VR* is a python function defining the volume ratio.  If it is not
     present, the volume ratio is 1.
+
+    *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the
+    C source code for the body of the volume, Iq, and Iqxy functions
+    respectively.  These can also be defined in the last source file.
 
 An *info* dictionary is constructed from the kernel meta data and
@@ -156,4 +174,6 @@
 __all__ = ["make, doc", "sources"]
 
+import sys
+import os
 import os.path
 
@@ -226,7 +246,7 @@
 #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
+// Standard mathematical constants:
+//   M_E, M_LOG2E, M_LOG10E, M_LN2, M_LN10, M_PI, M_PI_2=pi/2, M_PI_4=pi/4,
+//   M_1_PI=1/pi, M_2_PI=2/pi, M_2_SQRTPI=2/sqrt(pi), SQRT2, SQRT1_2=sqrt(1/2)
 // OpenCL defines M_constant_F for float constants, and nothing if double
 // is not enabled on the card, which is why these constants may be missing
@@ -256,4 +276,5 @@
     'qinit': "const real qi = q[i];",
     'qcall': "qi",
+    'qwork': ["q"],
     }
 
@@ -263,4 +284,5 @@
     'qinit': "const real qxi = qx[i];\n    const real qyi = qy[i];",
     'qcall': "qxi, qyi",
+    'qwork': ["qx", "qy"],
     }
 
@@ -318,5 +340,6 @@
 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];"""
+  const real %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\
+"""
 
 # Polydispersity loop body.
@@ -327,8 +350,20 @@
 const real weight = %(weight_product)s;
 if (weight > cutoff) {
-  ret += weight*%(fn)s(%(qcall)s, %(pcall)s);
-  norm += weight;
-  %(volume_norm)s
-}"""
+  const real I = %(fn)s(%(qcall)s, %(pcall)s);
+  if (I>=REAL(0.0)) { // scattering cannot be negative
+    ret += weight*I;
+    norm += weight;
+    %(volume_norm)s
+  }
+  //else { printf("exclude qx,qy,I:%%g,%%g,%%g\\n",%(qcall)s,I); }
+}
+//else { printf("exclude weight:%%g\\n",weight); }\
+"""
+
+# Use this when integrating over orientation
+SPHERICAL_CORRECTION="""\
+// Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta
+real spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*phi)) : REAL(1.0));\
+"""
 
 # Volume normalization.
@@ -337,6 +372,16 @@
 # a normalized weight.
 VOLUME_NORM="""const real vol_weight = %(weight)s;
-  vol += vol_weight*form_volume(%(pars)s);
-  norm_vol += vol_weight;"""
+    vol += vol_weight*form_volume(%(pars)s);
+    norm_vol += vol_weight;\
+"""
+
+# functions defined as strings in the .py module
+WORK_FUNCTION="""\
+real %(name)s(%(pars)s);
+real %(name)s(%(pars)s)
+{
+%(body)s
+}\
+"""
 
 # Documentation header for the module, giving the model name, its short
@@ -393,4 +438,5 @@
     vol_pars = info['partype']['volume']
     q_pars = KERNEL_2D if is_2D else KERNEL_1D
+    fn = q_pars['fn']
 
     # Build polydispersity loops
@@ -426,13 +472,19 @@
     fq_pars = [p[0] for p in info['parameters'][len(COMMON_PARAMETERS):]
                if p[0] in set(fixed_pars+pd_pars)]
+    if False and "phi" in pd_pars:
+        spherical_correction = [indent(SPHERICAL_CORRECTION, depth)]
+        weights = [p+"_w" for p in pd_pars]+['spherical_correction']
+    else:
+        spherical_correction = []
+        weights = [p+"_w" for p in pd_pars]
     subst = {
-        'weight_product': "*".join(p+"_w" for p in pd_pars),
+        'weight_product': "*".join(weights),
         'volume_norm': volume_norm,
-        'fn': q_pars['fn'],
+        'fn': fn,
         'qcall': q_pars['qcall'],
         'pcall': ", ".join(fq_pars), # skip scale and background
         }
     loop_body = [indent(LOOP_BODY%subst, depth)]
-    loops = "\n".join(loop_head+loop_body+loop_end)
+    loops = "\n".join(loop_head+spherical_correction+loop_body+loop_end)
 
     # declarations for non-pd followed by pd pars
@@ -465,4 +517,15 @@
         }
     kernel = KERNEL_TEMPLATE%subst
+
+    # If the working function is defined in the kernel metadata as a
+    # string, translate the string to an actual function definition
+    # and put it before the kernel.
+    if info[fn]:
+        subst = {
+            'name': fn,
+            'pars': ", ".join("real "+p for p in q_pars['qwork']+fq_pars),
+            'body': info[fn],
+            }
+        kernel = "\n".join((WORK_FUNCTION%subst, kernel))
     return kernel
 
@@ -514,5 +577,5 @@
     search_path = [ dirname(info['filename']),
                     abspath(joinpath(dirname(__file__),'models')) ]
-    return [_search(search_path) for f in info['source']]
+    return [_search(search_path, f) for f in info['source']]
 
 def make_model(info):
@@ -521,7 +584,18 @@
     found in the given search path.
     """
+    source = [open(f).read() for f in sources(info)]
+    # If the form volume is defined as a string, then wrap it in a
+    # function definition and place it after the external sources but
+    # before the kernel functions.  If the kernel functions are strings,
+    # they will be translated in the make_kernel call.
+    if info['form_volume']:
+        subst = {
+            'name': "form_volume",
+            'pars': ", ".join("real "+p for p in info['partype']['volume']),
+            'body': info['form_volume'],
+            }
+        source.append(WORK_FUNCTION%subst)
     kernel_Iq = make_kernel(info, is_2D=False)
     kernel_Iqxy = make_kernel(info, is_2D=True)
-    source = [open(f).read() for f in sources(info)]
     kernel = "\n\n".join([KERNEL_HEADER]+source+[kernel_Iq, kernel_Iqxy])
     return kernel
@@ -573,10 +647,12 @@
         name = kernel_module.name,
         title = kernel_module.title,
-        source = kernel_module.source,
         description = kernel_module.description,
         parameters = COMMON_PARAMETERS + kernel_module.parameters,
-        ER = getattr(kernel_module, 'ER', None),
-        VR = getattr(kernel_module, 'VR', None),
+        source = getattr(kernel_module, 'source', []),
         )
+    # Fill in attributes which default to None
+    info.update((k,getattr(kernel_module, k, None))
+                for k in ('ER', 'VR', 'form_volume', 'Iq', 'Iqxy'))
+    # Fill in the derived attributes
     info['limits'] = dict((p[0],p[3]) for p in info['parameters'])
     info['partype'] = categorize_parameters(info['parameters'])
@@ -596,65 +672,4 @@
     return DOC_HEADER%subst
 
-# Compiler platform details
-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(info):
-    """
-    Path to the compiled model defined by *info*.
-    """
-    from os.path import join as joinpath, split as splitpath, splitext
-    basename = splitext(splitpath(info['filename'])[1])[0]
-    return joinpath(DLL_PATH, basename+'.so')
-
-
-def dll_model(kernel_module, dtype=None):
-    """
-    Load the compiled model defined by *kernel_module*.
-
-    Recompile if any files are newer than the model file.
-
-    *dtype* is ignored.  Compiled files are always double.
-
-    The DLL is not loaded until the kernel is called so models an
-    be defined without using too many resources.
-    """
-    import tempfile
-    from sasmodels import dll
-
-    source, info = make(kernel_module)
-    source_files = sources(info) + [info['filename']]
-    newest = max(os.path.getmtime(f) for f in source_files)
-    dllpath = dll_path(info)
-    if not os.path.exists(dllpath) or os.path.getmtime(dllpath)<newest:
-        # Replace with a proper temp file
-        srcfile = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name'])
-        open(srcfile, 'w').write(source)
-        os.system(COMPILE%(srcfile, dllpath))
-        ## comment the following to keep the generated c file
-        os.unlink(srcfile)
-    return dll.DllModel(dllpath, info)
-
-
-def opencl_model(kernel_module, dtype="single"):
-    """
-    Load the OpenCL model defined by *kernel_module*.
-
-    Access to the OpenCL device is delayed until the kernel is called
-    so models can be defined without using too many resources.
-    """
-    from sasmodels import gpu
-
-    source, info = make(kernel_module)
-    ## 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, info, dtype)
 
 
Index: sasmodels/gpu.py
===================================================================
--- sasmodels/gpu.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/gpu.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -25,6 +25,4 @@
 
 """
-import warnings
-
 import numpy as np
 import pyopencl as cl
@@ -32,6 +30,4 @@
 
 from . import gen
-
-from .gen import F32, F64
 
 F32_DEFS = """\
@@ -52,5 +48,18 @@
 # larger than necessary given that cost grows as npts^k where k is the number
 # of polydisperse parameters.
-MAX_LOOPS = 1024
+MAX_LOOPS = 2048
+
+def load_model(kernel_module, dtype="single"):
+    """
+    Load the OpenCL model defined by *kernel_module*.
+
+    Access to the OpenCL device is delayed until the kernel is called
+    so models can be defined without using too many resources.
+    """
+    source, info = gen.make(kernel_module)
+    ## for debugging, save source to a .cl file, edit it, and reload as model
+    open(info['name']+'.cl','w').write(source)
+    #source = open(info['name']+'.cl','r').read()
+    return GpuModel(source, info, dtype)
 
 ENV = None
@@ -103,8 +112,8 @@
     """
     dtype = np.dtype(dtype)
-    if dtype==F64 and not all(has_double(d) for d in context.devices):
+    if dtype==gen.F64 and not all(has_double(d) for d in context.devices):
         raise RuntimeError("Double precision not supported for devices")
 
-    header = F64_DEFS if dtype == F64 else F32_DEFS
+    header = F64_DEFS if dtype == gen.F64 else F32_DEFS
     # Note: USE_SINCOS makes the intel cpu slower under opencl
     if context.devices[0].type == cl.device_type.GPU:
@@ -158,5 +167,5 @@
     is an optional extension which may not be available on all devices.
     """
-    def __init__(self, source, info, dtype=F32):
+    def __init__(self, source, info, dtype=gen.F32):
         self.info = info
         self.source = source
@@ -221,5 +230,5 @@
     buffer will be released when the data object is freed.
     """
-    def __init__(self, q_vectors, dtype=F32):
+    def __init__(self, q_vectors, dtype=gen.F32):
         env = environment()
         self.nq = q_vectors[0].size
@@ -273,5 +282,5 @@
         env = environment()
         self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                  MAX_LOOPS*input.dtype.itemsize)
+                                  2*MAX_LOOPS*input.dtype.itemsize)
                         for _ in env.queues]
         self.res_b = [cl.Buffer(env.context, mf.READ_WRITE,
@@ -281,14 +290,15 @@
 
     def __call__(self, pars, pd_pars, cutoff=1e-5):
-        real = np.float32 if self.input.dtype == F32 else np.float64
+        real = np.float32 if self.input.dtype == gen.F32 else np.float64
         fixed = [real(p) for p in pars]
         cutoff = real(cutoff)
         loops = np.hstack(pd_pars)
         loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten()
-        loops_N = [np.uint32(len(p[0])) for p in pd_pars]
+        Nloops = [np.uint32(len(p[0])) for p in pd_pars]
+        #print "loops",Nloops, loops
 
         #import sys; print >>sys.stderr,"opencl eval",pars
         #print "opencl eval",pars
-        if len(loops) > MAX_LOOPS:
+        if len(loops) > 2*MAX_LOOPS:
             raise ValueError("too many polydispersity points")
         device_num = 0
@@ -300,5 +310,5 @@
         #ctx = environment().context
         #loops_bi = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops)
-        args = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + loops_N
+        args = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + Nloops
         self.kernel(queuei, self.input.global_size, None, *args)
         cl.enqueue_copy(queuei, self.res, res_bi)
Index: sasmodels/models/README.rst
===================================================================
--- sasmodels/models/README.rst	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/models/README.rst	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -7,23 +7,12 @@
 model.
 
-cylindermodel_onefile.py is cylinder.c+cylinder.py merged into one file.
-This doesn't actually run yet since sasmodels.gen has not been updated
-to support it.  It exists as a proposal.  Note that the function declaration
-has been removed since there is enough information in the parameter
-definitions to generate it automatically.  Note also that "source" which
-used to be all the source has been renamed "includes".
-
-One-file models could coexist with the py+c file models by checking for the
-existence of c_blah and creating the appropriate function wrappers.  These
-would be appended after any include files.  You shouldn't mix the two forms
-within a single model since form_volume needs to be defined before
-Iq/Iqxy but after the libraries.
+lamellar.py is an example of a single file model with embedded C code.
 
 Note: may want to rename form_volume to calc_volume and Iq/Iqxy to
-calc_Iq/calc_Iqxy in both the c+py and the one file forms so that the
-names are more predictable.  Similarly ER/VR go to calc_ER/calc_VR.
+calc_Iq/calc_Iqxy. Similarly ER/VR go to calc_ER/calc_VR.
 
 Note: It is possible to translate python code automatically to opencl, using
-something like numba, clyther, shedskin or pypy.
+something like numba, clyther, shedskin or pypy, so maybe the kernel functions
+could be implemented without any C syntax.
 
 Magnetism hasn't been implemented yet.  We may want a separate Imagnetic
@@ -35,3 +24,4 @@
 
 Need to write code to generate the polydispersity loops in python for
-kernels that are only implemented in python.
+kernels that are only implemented in python.  Also, need to implement
+an example kernel directly in python.
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/capped_cylinder.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,149 @@
+real form_volume(real radius, real cap_radius, real length);
+real Iq(real q, real sld, real solvent_sld,
+    real radius, real cap_radius, real length);
+real Iqxy(real qx, real qy, real sld, real solvent_sld,
+    real radius, real cap_radius, real length, real theta, real phi);
+
+// Integral over a convex lens kernel for t in [h/R,1].  See the docs for
+// the definition of the function being integrated.
+//   q is the magnitude of the q vector.
+//   h is the length of the lens "inside" the cylinder.  This negative wrt the
+//       definition of h in the docs.
+//   cap_radius is the radius of the lens
+//   length is the cylinder length, or the separation between the lens halves
+//   alpha is the angle of the cylinder wrt q.
+real _cap_kernel(real q, real h, real cap_radius, real length,
+                 real sin_alpha, real cos_alpha);
+real _cap_kernel(real q, real h, real cap_radius, real length,
+                 real sin_alpha, real cos_alpha)
+{
+    // For speed, we are pre-calculating terms which are constant over the
+    // kernel.
+    const real upper = REAL(1.0);
+    const real lower = h/cap_radius; // integral lower bound
+    const real m = q*cos_alpha*cap_radius; // cos argument slope
+    const real b = q*cos_alpha*(REAL(0.5)*length-h); // cos argument intercept
+    const real qrst = q*sin_alpha*cap_radius; // Q*R*sin(theta)
+    real total = REAL(0.0);
+    for (int i=0; i<76 ;i++) {
+        // translate a point in [-1,1] to a point in [lower,upper]
+        //const real t = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        const real t = REAL(0.5)*(Gauss76Z[i]*(upper-lower)+upper+lower);
+        const real radical = REAL(1.0) - t*t;
+        const real caparg = qrst*sqrt(radical); // cap bessel function arg
+        const real be = (caparg == REAL(0.0) ? REAL(0.5) : J1(caparg)/caparg);
+        const real Fq = cos(m*t + b) * radical * be;
+        total += Gauss76Wt[i] * Fq;
+    }
+    // translate dx in [-1,1] to dx in [lower,upper]
+    //const real form = (upper-lower)/2.0*total;
+    const real integral = REAL(0.5)*(upper-lower)*total;
+    return REAL(4.0)*M_PI*cap_radius*cap_radius*cap_radius*integral;
+}
+
+real form_volume(real radius, real cap_radius, real length)
+{
+    // cap radius should never be less than radius when this is called
+    // Note: cap volume = pi hc/6 * (3 a^2 + hc^2), where a is the cylinder
+    // radius and hc is the height of the cap.  Multiply by two for both ends.
+    // So:
+    //      cap V = pi hc (r^2 + hc^2/3)
+    //      cylinder V = pi r^2 L
+    //      V = cylinder V + cap V
+    //        = pi r^2 L + pi hc (r^2 + hc^2/3)
+    //        = pi * (r^2 (L+hc) + hc^3/3)
+    const real hc = cap_radius - sqrt(cap_radius*cap_radius - radius*radius);
+    return M_PI*(radius*radius*(length+hc) + REAL(0.333333333333333)*hc*hc*hc);
+}
+
+real Iq(real q,
+    real sld,
+    real solvent_sld,
+    real radius,
+    real cap_radius,
+    real length)
+{
+    real sn, cn; // slots to hold sincos function output
+
+    // Exclude invalid inputs.
+    if (cap_radius < radius) return REAL(-1.0);
+
+    const real lower = REAL(0.0);
+    const real upper = M_PI_2;
+    const real h = sqrt(cap_radius*cap_radius - radius*radius); // negative h
+    real total = REAL(0.0);
+    for (int i=0; i<76 ;i++) {
+        // translate a point in [-1,1] to a point in [lower,upper]
+        const real alpha= REAL(0.5)*(Gauss76Z[i]*(upper-lower) + upper + lower);
+        SINCOS(alpha, sn, cn);
+
+        const real cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn);
+
+        // 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 multiply and divide sin(alpha).
+        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 cyl_Fq = M_PI*radius*radius*length*REAL(2.0)*bj*si;
+
+        // Volume weighted average F(q)
+        const real Aq = cyl_Fq + cap_Fq;
+        total += Gauss76Wt[i] * Aq * Aq * sn; // sn for spherical coord integration
+    }
+    // translate dx in [-1,1] to dx in [lower,upper]
+    const real form = total * REAL(0.5)*(upper-lower);
+
+    // 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);
+    return REAL(1.0e-4) * form * s * s; // form_volume(radius, cap_radius, length);
+}
+
+
+real Iqxy(real qx, real qy,
+    real sld,
+    real solvent_sld,
+    real radius,
+    real cap_radius,
+    real length,
+    real theta,
+    real phi)
+{
+    real sn, cn; // slots to hold sincos function output
+
+    // Exclude invalid inputs.
+    if (cap_radius < radius) return REAL(-1.0);
+
+    // 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 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); // rod angle relative to q
+    SINCOS(alpha, sn, cn);
+
+    const real h = sqrt(cap_radius*cap_radius - radius*radius); // negative h
+    const real cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn);
+
+    // 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 multiply and divide sin(alpha).
+    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 cyl_Fq = M_PI*radius*radius*length*REAL(2.0)*bj*si;
+
+    // Volume weighted average F(q)
+    const real Aq = cyl_Fq + cap_Fq;
+
+    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
+    const real s = (sld - solvent_sld);
+    return REAL(1.0e-4) * Aq * Aq * s * s; // form_volume(radius, cap_radius, length);
+}
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/capped_cylinder.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,143 @@
+r"""
+Calculates the scattering from a cylinder with spherical section end-caps.
+This model simply becomes the `convex_lens`_ when the length of the cylinder
+$L=0$, that is, a sphereocylinder with end caps that have a radius larger
+than that of the cylinder and the center of the end cap radius lies within
+the cylinder. See the diagram for the details of the geometry and
+restrictions on parameter values.
+
+Definitions
+-----------
+
+The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale.
+
+The capped cylinder geometry is defined as
+
+.. image:: img/image112.JPG
+
+where $r$ is the radius of the cylinder. All other parameters are as defined
+in the diagram. Since the end cap radius $R \ge r$ and by definition for this
+geometry $h < 0$, $h$ is then defined by $r$ and $R$ as
+
+.. math::
+
+    h = - \sqrt{R^2 - r^2}
+
+The scattered intensity $I(q)$ is calculated as
+
+.. math::
+
+    I(q) = \frac{(\Delta \rho)^2}{V} left<A^2(q)\right>
+
+where the amplitude $A(q)$ is given as
+
+.. math::
+
+    A(Q) = \pi r^2L \frac{sin(Q(L/2) \cos \theta)}{Q(L/2) \cos \theta} \
+        \frac{2 J_1(Q r \sin \theta)}{Q r \sin \theta} \
+        + 4 \pi R^3 \int_{-h/R}^1{dt \cos [ Q cos\theta (Rt + h + L/2)] \
+        x (1-t^2)\frac{J_1[Q R \sin \theta (1-t^2)^{1/2}]}{Q R \sin \theta (1-t^2)^{1/2}
+
+The $\left< \ldots \right>$ brackets denote an average of the structure over
+all orientations. $\left< A^2(q)\right>$ is then the form factor, $P(q)$.
+The scale factor is equivalent to the volume fraction of cylinders, each of
+volume, $V$. Contrast is the difference of scattering length densities of
+the cylinder and the surrounding solvent.
+
+The volume of the capped cylinder is (with $h$ as a positive value here)
+
+.. math::
+
+    V = \pi r_c^2 L + \frac{2\pi}{3}(R-h)^2(2R + h)
+
+
+and its radius-of-gyration is
+
+    R_g^2 = \left[ \frac{12}{5}R^5 + R^4(6h+\frac{3}{2}L) \
+        + R^2(4h^2 + L^2 + 4Lh) + R^2(3Lh^2 + [frac{3}{2}L^2h) \
+        + \frac{2}{5}h^5 - \frac{1}{2}Lh^4 - \frac{1}{2}L^2h^3 \
+        + \frac{1}{4}L^3r^2 + \frac{3}{2}Lr^4\right] \
+        (4R^3 6R^2h - 2h^3 + 3r^2L)^{-1}
+
+
+**The requirement that $R \ge r$ is not enforced in the model! It is up to
+you to restrict this during analysis.**
+
+Figure :num:`figure #capped-cylinder-1d` shows the output produced by
+a running the 1D capped cylinder model, using *qmin* = 0.001 |Ang^-1|,
+*qmax* = 0.7 |Ang^-1| and  the default values of the parameters.
+
+.. _capped-cylinder-1d:
+
+.. figure:: img/image117.jpg
+
+    1D plot using the default values (w/256 data point).
+
+The 2D scattering intensity is calculated similar to the 2D cylinder model.
+Figure :num:`figure #capped-cylinder-2d` shows the output for $\theta=45^\circ$
+and $\phi=0^\circ$ with default values for the other parameters.
+
+.. _capped-cylinder-2d:
+
+.. figure:: img/image118.JPG
+
+    2D plot (w/(256X265) data points).
+
+.. figure:: img/image061.JPG
+
+    Definition of the angles for oriented 2D cylinders.
+
+.. figure:: img/image062.jpg
+
+    Examples of the angles for oriented pp against the detector plane.
+
+REFERENCE
+
+H Kaya, *J. Appl. Cryst.*, 37 (2004) 223-230
+
+H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda and errata)
+"""
+from numpy import inf
+
+name = "capped_cylinder"
+title = "Right circular cylinder with spherical end caps and uniform SLD"
+description = """That is, a sphereocylinder
+	with end caps that have a radius larger than
+	that of the cylinder and the center of the
+	end cap radius lies within the cylinder.
+	Note: As the length of cylinder -->0,
+	it becomes a ConvexLens.
+	It must be that radius <(=) cap_radius.
+	[Parameters];
+	scale: volume fraction of spheres,
+	background:incoherent background,
+	radius: radius of the cylinder,
+	length: length of the cylinder,
+	cap_radius: radius of the semi-spherical cap,
+	sld: SLD of the capped cylinder,
+	solvent_sld: SLD of the solvent.
+"""
+
+parameters = [
+    #   [ "name", "units", default, [lower, upper], "type",
+    #     "description" ],
+    [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "",
+      "Cylinder scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "",
+      "Solvent scattering length density" ],
+    [ "radius", "Ang",  20, [0, inf], "volume",
+      "Cylinder radius" ],
+    [ "cap_radius", "Ang", 20, [0, inf], "volume",
+      "Cap radius" ],
+    [ "length", "Ang",  400, [0, inf], "volume",
+      "Cylinder length" ],
+    [ "theta", "degrees", 60, [-inf, inf], "orientation",
+      "In plane angle" ],
+    [ "phi", "degrees", 60, [-inf, inf], "orientation",
+      "Out of plane angle" ],
+    ]
+
+source = [ "lib/J1.c", "lib/gauss76.c", "capped_cylinder.c" ]
+
+
+
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/core_shell_cylinder.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,92 @@
+real form_volume(real radius, real thickness, real length);
+real Iq(real q, real core_sld, real shell_sld, real solvent_sld,
+    real radius, real thickness, real length);
+real Iqxy(real qx, real qy, real core_sld, real shell_sld, real solvent_sld,
+    real radius, real thickness, real length, real theta, real phi);
+
+// twovd = 2 * volume * delta_rho
+// besarg = q * R * sin(alpha)
+// siarg = q * L/2 * cos(alpha)
+real _cyl(real twovd, real besarg, real siarg);
+real _cyl(real twovd, real besarg, real siarg)
+{
+    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 twovd*si*bj;
+}
+
+real form_volume(real radius, real thickness, real length)
+{
+    return M_PI*(radius+thickness)*(radius+thickness)*(length+2*thickness);
+}
+
+real Iq(real q,
+    real core_sld,
+    real shell_sld,
+    real solvent_sld,
+    real radius,
+    real thickness,
+    real length)
+{
+    // precalculate constants
+    const real core_qr = q*radius;
+    const real core_qh = q*REAL(0.5)*length;
+    const real core_twovd = REAL(2.0) * form_volume(radius,0,length)
+                            * (core_sld-shell_sld);
+    const real shell_qr = q*(radius + thickness);
+    const real shell_qh = q*(REAL(0.5)*length + thickness);
+    const real shell_twovd = REAL(2.0) * form_volume(radius,thickness,length)
+                             * (shell_sld-solvent_sld);
+    real total = REAL(0.0);
+    // real lower=0, upper=M_PI_2;
+    for (int i=0; i<76 ;i++) {
+        // translate a point in [-1,1] to a point in [lower,upper]
+        //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        real sn, cn;
+        const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);
+        SINCOS(alpha, sn, cn);
+        const real fq = _cyl(core_twovd, core_qr*sn, core_qh*cn)
+            + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn);
+        total += Gauss76Wt[i] * fq * fq * sn;
+    }
+    // translate dx in [-1,1] to dx in [lower,upper]
+    //const real form = (upper-lower)/2.0*total;
+    return REAL(1.0e-4) * total * M_PI_4;
+}
+
+
+real Iqxy(real qx, real qy,
+    real core_sld,
+    real shell_sld,
+    real solvent_sld,
+    real radius,
+    real thickness,
+    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);
+
+    const real core_qr = q*radius;
+    const real core_qh = q*REAL(0.5)*length;
+    const real core_twovd = REAL(2.0) * form_volume(radius,0,length)
+                            * (core_sld-shell_sld);
+    const real shell_qr = q*(radius + thickness);
+    const real shell_qh = q*(REAL(0.5)*length + thickness);
+    const real shell_twovd = REAL(2.0) * form_volume(radius,thickness,length)
+                             * (shell_sld-solvent_sld);
+
+    SINCOS(alpha, sn, cn);
+    const real fq = _cyl(core_twovd, core_qr*sn, core_qh*cn)
+        + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn);
+    return REAL(1.0e-4) * fq * fq;
+}
Index: sasmodels/models/core_shell_cylinder.py
===================================================================
--- sasmodels/models/core_shell_cylinder.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/core_shell_cylinder.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,164 @@
+# core shell cylinder model
+# Note: model title and parameter table are inserted automatically
+r"""
+The form factor is normalized by the particle volume.
+
+Definition
+----------
+
+The output of the 2D scattering intensity function for oriented core-shell
+cylinders is given by (Kline, 2006)
+
+.. math::
+
+    P(q,\alpha) = \frac{\text{scale}}{V_s} f^2(q) + \text{background}
+
+where
+
+.. math::
+
+    f(q) = (\rho_c - \rho_s) V_c \
+           \frac{\sin( Q L/2 \cos \alpha)}{Q L/2 \cos \alpha} \
+           \frac{2 J_1 (Q R \sin \alpha)}{Q R \sin \alpha}
+           + (\rho_s - \rho_\text{solv}) V_s \
+           \frac{\sin( Q (L/2+T) \cos \alpha)}{Q (L/2+T) \cos \alpha}
+           \frac{2 J_1 (Q (R+T) \sin \alpha)}{Q (R+T) \sin \alpha}
+
+and
+
+.. math::
+
+    V_s = \pi (R + T)^2 \dot (L + 2T)
+
+
+and $\alpha$ is the angle between the axis of the cylinder and $\vec q$,
+$V_s$ is the volume of the outer shell (i.e. the total volume, including
+the shell), $V_c$ is the volume of the core, $L$ is the length of the core,
+$R$ is the radius of the core, $T$ is the thickness of the shell, $\rho_c$
+is the scattering length density of the core, $\rho_s$ is the scattering
+length density of the shell, $\rho_\text{solv}$ is the scattering length
+density of the solvent, and *background* is the background level.  The outer
+radius of the shell is given by $R+T$ and the total length of the outer
+shell is given by $L+2T$. $J1$ is the first order Bessel function.
+
+.. figure:: img/image069.JPG
+
+    Core shell cylinder schematic.
+
+To provide easy access to the orientation of the core-shell cylinder, we
+define the axis of the cylinder using two angles $\theta$ and $\phi$. As
+for the case of the cylinder, those angles are defined in
+Figure :num:`figure #cylinder-orientation`.
+
+NB: The 2nd virial coefficient of the cylinder is calculated based on
+the radius and 2 length values, and used as the effective radius for
+$S(Q)$ when $P(Q) \dot S(Q)$ is applied.
+
+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
+----------
+
+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 #core-shell-cylinder-comparison-1d` shows a comparison
+of the 1D output of our model and the output of the NIST software.
+
+.. _core-shell-cylinder-comparison-1d:
+
+.. figure:: img/image070.JPG
+
+    Comparison of the SasView scattering intensity for a core-shell cylinder
+    with the output of the NIST SANS analysis software. The parameters were
+    set to: *scale=1.0 |Ang|*, *radius=20 |Ang|*, *thickness=10 |Ang|*,
+    *length=400 |Ang|*, *core_sld=1e-6 |Ang^-2|*, *shell_sld=4e-6 |Ang^-2|*,
+    *solvent_sld=1e-6 |Ang^-2|*, and *background=0.01 |cm^-1|*.
+
+Averaging over a distribution of orientation is done by evaluating the
+equation above. 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 #core-shell-cylinder-comparison-2d` shows the result
+of such a cross-check.
+
+.. _core-shell-cylinder-comparison-2d:
+
+.. figure:: img/image071.JPG
+
+    Comparison of the intensity for uniformly distributed core-shell
+    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|*, *thickness=10 |Ang|*, *length*=400 |Ang|*,
+    *core_sld=1e-6 |Ang^-2|*, *shell_sld=4e-6 |Ang^-2|*,
+    *solvent_sld=1e-6 |Ang^-2|*, and *background=0.0 |cm^-1|*.
+
+2013/11/26 - Description reviewed by Heenan, R.
+"""
+
+from numpy import pi, inf
+
+name = "cylinder"
+title = "Right circular cylinder with a core-shell scattering length density profile."
+description = """
+P(q,alpha)= scale/Vs*f(q)^(2) + background,
+      where: f(q)= 2(core_sld - solvant_sld)
+		* Vc*sin[qLcos(alpha/2)]
+		/[qLcos(alpha/2)]*J1(qRsin(alpha))
+		/[qRsin(alpha)]+2(shell_sld-solvent_sld)
+		*Vs*sin[q(L+T)cos(alpha/2)][[q(L+T)
+		*cos(alpha/2)]*J1(q(R+T)sin(alpha))
+		/q(R+T)sin(alpha)]
+
+	alpha:is the angle between the axis of
+		the cylinder and the q-vector
+	Vs: the volume of the outer shell
+	Vc: the volume of the core
+	L: the length of the core
+    	shell_sld: the scattering length density of the shell
+	solvent_sld: the scattering length density of the solvent
+	background: the background
+	T: the thickness
+    	R+T: is the outer radius
+ 	L+2T: The total length of the outershell
+	J1: the first order Bessel function
+ 	theta: axis_theta of the cylinder
+ 	phi: the axis_phi of the cylinder
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "core_sld", "1e-6/Ang^2", 4, [-inf,inf], "",
+      "Cylinder core scattering length density" ],
+    [ "shell_sld", "1e-6/Ang^2", 4, [-inf,inf], "",
+      "Cylinder shell scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "",
+      "Solvent scattering length density" ],
+    [ "radius", "Ang",  20, [0, inf], "volume",
+      "Cylinder core radius" ],
+    [ "thickness", "Ang",  20, [0, inf], "volume",
+      "Cylinder shell thickness" ],
+    [ "length", "Ang",  400, [0, inf], "volume",
+      "Cylinder length" ],
+    [ "theta", "degrees", 60, [-inf, inf], "orientation",
+      "In plane angle" ],
+    [ "phi", "degrees", 60, [-inf, inf], "orientation",
+      "Out of plane angle" ],
+    ]
+
+source = [ "lib/J1.c", "lib/gauss76.c", "core_shell_cylinder.c"]
+
+def ER(radius, thickness, length):
+    radius = radius + thickness
+    length = length + 2*thickness
+    ddd = 0.75*radius*(2*radius*length + (length+radius)*(length+pi*radius))
+    return 0.5 * (ddd)**(1./3.)
+
+def VR(radius, thickness, length):
+    whole = pi * (radius+thickness)**2 * (length+2*thickness)
+    core = pi * radius**2 * length
+    return whole, whole-core
+
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/models/cylinder.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -1,5 +1,17 @@
 real form_volume(real radius, real length);
 real Iq(real q, real sld, real solvent_sld, real radius, real length);
-real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi);
+real Iqxy(real qx, real qy, real sld, real solvent_sld,
+    real radius, real length, real theta, real phi);
+
+// twovd = 2 * volume * delta_rho
+// besarg = q * R * sin(alpha)
+// siarg = q * L/2 * cos(alpha)
+real _cyl(real twovd, real besarg, real siarg);
+real _cyl(real twovd, real besarg, real siarg)
+{
+    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 twovd*si*bj;
+}
 
 real form_volume(real radius, real length)
@@ -14,22 +26,21 @@
     real length)
 {
-    const real halflength = REAL(0.5)*length;
-    real summ = REAL(0.0);
+    const real qr = q*radius;
+    const real qh = q*REAL(0.5)*length;
+    const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length);
+    real total = REAL(0.0);
     // real lower=0, upper=M_PI_2;
     for (int i=0; i<76 ;i++) {
         // translate a point in [-1,1] to a point in [lower,upper]
-        //const real zi = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
-        const real zi = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);
-        summ += Gauss76Wt[i] * CylKernel(q, radius, halflength, zi);
+        //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);
+        real sn, cn;
+        SINCOS(alpha, sn, cn);
+        const real fq = _cyl(twovd, qr*sn, qh*cn);
+        total += Gauss76Wt[i] * fq * fq * sn;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
-    //const real form = (upper-lower)/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) * form_volume(radius, length);
-    return REAL(1.0e-4) * form * s * s;
+    //const real form = (upper-lower)/2.0*total;
+    return REAL(1.0e-4) * total * M_PI_4;
 }
 
@@ -43,4 +54,7 @@
     real phi)
 {
+    // TODO: check that radius<0 and length<0 give zero scattering.
+    // This should be the case since the polydispersity weight vector should
+    // be zero length, and this function never called.
     real sn, cn; // slots to hold sincos function output
 
@@ -54,19 +68,7 @@
     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).
+    const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length);
     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) * form_volume(radius, length);
-    return REAL(1.0e-4) * form * s * s; // * correction;
+    const real fq = _cyl(twovd, q*radius*sn, q*REAL(0.5)*length*cn);
+    return REAL(1.0e-4) * fq * fq; // * correction;
 }
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision a7684e5745716d713f91b1ebc1027632a3ff9d36)
+++ sasmodels/models/cylinder.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -1,2 +1,3 @@
+# cylinder model
 # Note: model title and parameter table are inserted automatically
 r"""
@@ -31,7 +32,7 @@
 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`.
+are defined in Figure :num:`figure #cylinder-orientation`.
 
-.. _CylinderModel-orientation:
+.. _cylinder-orientation:
 
 .. figure:: img/image061.JPG   (should be img/cylinder-1.jpg, or img/cylinder-orientation.jpg)
@@ -64,8 +65,8 @@
 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
+Figure :num:`figure #cylinder-compare` shows a comparison of
 the 1D output of our model and the output of the NIST software.
 
-.. _CylinderModel-compare:
+.. _cylinder-compare:
 
 .. figure:: img/image065.JPG
@@ -91,8 +92,8 @@
 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
+Figure :num:`figure #cylinder-crosscheck` shows the result of
 such a cross-check.
 
-.. _CylinderModel-crosscheck:
+.. _cylinder-crosscheck:
 
 .. figure:: img/image066.JPG
@@ -111,5 +112,5 @@
 title = "Right circular cylinder with uniform scattering length density."
 description = """
-     f(q)= 2*(sldCyl - sldSolv)*V*sin(qLcos(alpha/2))
+     P(q)= 2*(sld - solvent_sld)*V*sin(qLcos(alpha/2))
             /[qLcos(alpha/2)]*J1(qRsin(alpha/2))/[qRsin(alpha)]
 
@@ -119,10 +120,10 @@
             L: Length of the cylinder
             J1: The bessel function
-            alpha: angle betweenthe axis of the
+            alpha: angle between the axis of the
             cylinder and the q-vector for 1D
             :the ouput is P(q)=scale/V*integral
             from pi/2 to zero of...
-            f(q)^(2)*sin(alpha)*dalpha+ bkg
-    """
+            f(q)^(2)*sin(alpha)*dalpha + background
+"""
 
 parameters = [
@@ -143,5 +144,5 @@
     ]
 
-source = [ "lib/J1.c", "lib/gauss76.c", "lib/cylkernel.c", "cylinder.c"]
+source = [ "lib/J1.c", "lib/gauss76.c", "cylinder.c" ]
 
 def ER(radius, length):
Index: sasmodels/models/cylinder_clone.c
===================================================================
--- sasmodels/models/cylinder_clone.c	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/models/cylinder_clone.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -2,4 +2,16 @@
 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);
+
+
+// twovd = 2 * volume * delta_rho
+// besarg = q * R * sin(alpha)
+// siarg = q * L/2 * cos(alpha)
+real _cyl(real twovd, real besarg, real siarg, real alpha);
+real _cyl(real twovd, real besarg, real siarg, real alpha)
+{
+    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 twovd*si*bj;
+}
 
 real form_volume(real radius, real length)
@@ -7,5 +19,4 @@
     return M_PI*radius*radius*length;
 }
-
 real Iq(real q,
     real sldCyl,
@@ -14,19 +25,21 @@
     real length)
 {
-    const real h = REAL(0.5)*length;
-    real summ = REAL(0.0);
+    const real qr = q*radius;
+    const real qh = q*REAL(0.5)*length;
+    const real twovd = REAL(2.0)*(sldCyl-sldSolv)*form_volume(radius, length);
+    real total = REAL(0.0);
+    // real lower=0, upper=M_PI_2;
     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);
+        // translate a point in [-1,1] to a point in [lower,upper]
+        //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);
+        real sn, cn;
+        SINCOS(alpha, sn, cn);
+        const real fq = _cyl(twovd, qr*sn, qh*cn, alpha);
+        total += Gauss76Wt[i] * fq * fq * sn;
     }
-    //const real form = (uplim-lolim)/2.0*summ;
-    const real form = summ * M_PI_4;
-
-    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
-    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
-    // The additional volume factor is for polydisperse volume normalization.
-    const real s = (sldCyl - sldSolv) * form_volume(radius, length);
-    return REAL(1.0e8) * form * s * s;
+    // translate dx in [-1,1] to dx in [lower,upper]
+    //const real form = (upper-lower)/2.0*total;
+    return REAL(1.0e8) * total * M_PI_4;
 }
 
@@ -50,19 +63,9 @@
     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).
+    const real qr = q*radius;
+    const real qh = q*REAL(0.5)*length;
+    const real twovd = REAL(2.0)*(sldCyl-sldSolv)*form_volume(radius, length);
     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 = (sldCyl - sldSolv) * form_volume(radius, length);
-    return REAL(1.0e8) * form * s * s * spherical_integration;
+    const real fq = _cyl(twovd, qr*sn, qh*cn, alpha);
+    return REAL(1.0e8) * fq * fq * spherical_integration;
 }
Index: sasmodels/models/cylinder_clone.py
===================================================================
--- sasmodels/models/cylinder_clone.py	(revision a7684e5745716d713f91b1ebc1027632a3ff9d36)
+++ sasmodels/models/cylinder_clone.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -150,5 +150,5 @@
       "Out of plane angle" ],
     ]
-source = [ "lib/J1.c", "lib/gauss76.c", "lib/cylkernel.c", "cylinder_clone.c"]
+source = [ "lib/J1.c", "lib/gauss76.c", "cylinder_clone.c" ]
 
 def ER(radius, length):
Index: sasmodels/models/cylinder_onefile.py
===================================================================
--- sasmodels/models/cylinder_onefile.py	(revision a7684e5745716d713f91b1ebc1027632a3ff9d36)
+++ sasmodels/models/cylinder_onefile.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -154,9 +154,9 @@
 source = [ "lib/J1.c", "lib/gauss76.c", "lib/cylkernel.c" ]
 
-c_form_volume = """
+form_volume = """
     return M_PI*radius*radius*length;
-"""
-
-c_Iq = """
+    """
+
+Iq = """
     const real h = REAL(0.5)*length;
     real summ = REAL(0.0);
@@ -174,7 +174,7 @@
     const real s = (sld - solvent_sld) * form_volume(radius, length);
     return REAL(1.0e-4) * form * s * s;
-"""
-
-c_Iqxy = """
+    """
+
+Iqxy = """
     real sn, cn; // slots to hold sincos function output
 
@@ -204,5 +204,5 @@
     const real s = (sld - solvent_sld) * form_volume(radius, length);
     return REAL(1.0e-4) * form * s * s; // * correction;
-"""
+    """
 
 def ER(radius, length):
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
+++ sasmodels/models/ellipsoid.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -1,48 +1,22 @@
-/* PARAMETERS
+real form_volume(real rpolar, real requatorial);
+real Iq(real q, real sld, real solvent_sld, real rpolar, real requatorial);
+real Iqxy(real qx, real qy, real sld, real solvent_sld,
+    real rpolar, real requatorial, real theta, real phi);
+
+real _ellipsoid_kernel(real q, real rpolar, real requatorial, real cos_alpha);
+real _ellipsoid_kernel(real q, real rpolar, real requatorial, real cos_alpha)
 {
-name: "ellipsoid",
-title: "Ellipsoid with uniform scattering length density",
-include: [ "lib/gauss76.c" ],
-parameters: [
-   // [ "name", "units", default, [lower, upper], "type", "description" ],
-   [ "sld", "1e-6/Ang^2", 4, [-Infinity,Infinity], "",
-     "Cylinder scattering length density" ],
-   [ "solvent_sld", "1e-6/Ang^2", 1, [-Infinity,Infinity], "",
-     "Solvent scattering length density" ],
-   [ "a", "Ang",  20, [0, Infinity], "volume",
-     "Cylinder radius" ],
-   [ "b", "Ang",  20, [0, Infinity], "volume",
-     "Cylinder length" ],
-   [ "theta", "degrees", 60, [-Infinity, Infinity], "orientation",
-     "In plane angle" ],
-   [ "phi", "degrees", 60, [-Infinity, Infinity], "orientation",
-     "Out of plane angle" ],
-],
-}
-PARAMETERS END
-
-DOCUMENTATION
-.. _EllipseModel:
-
-DOCUMENTATION END
-*/
-
-real form_volume(real a, real b);
-real Iq(real qx, real qy, real sld, real solvent_sld, real a, real b);
-real Iqxy(real qx, real qy, real sld, real solvent_sld, real a, real b, real theta, real phi);
-
-real form_volume(real a, real b)
-{
-    return REAL(1.333333333333333)*M_PI_2*a*b*b;
+    real sn, cn;
+    real ratio = rpolar/requatorial;
+    const real u = q*requatorial*sqrt(REAL(1.0)
+                   + cos_alpha*cos_alpha*(ratio*ratio - REAL(1.0)));
+    SINCOS(u, sn, cn);
+    const real f = ( u==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-u*cn)/(u*u*u) );
+    return f*f;
 }
 
-real ellipsoid_kernel(double q, double b, double a, double dum)
+real form_volume(real rpolar, real requatorial)
 {
-    real sn, cn;
-    const real nu = a/b;
-    const real arg = q * b * sqrt(REAL(1.0)+(dum*dum*(nu*nu--REAL(1.0))));
-    SINCOS(arg, sn, cn);
-    const real f = (arg==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-arg*cn)/(arg*arg*arg);
-    return f*f;
+    return REAL(1.333333333333333)*M_PI*rpolar*requatorial*requatorial;
 }
 
@@ -50,16 +24,18 @@
     real sld,
     real solvent_sld,
-    real a,
-    real b)
+    real rpolar,
+    real requatorial)
 {
-    real summ = REAL(0.0);
+    //const real lower = REAL(0.0);
+    //const real upper = REAL(1.0);
+    real total = REAL(0.0);
     for (int i=0;i<76;i++) {
-        //const real zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
-        zi = ( Gauss76Z[i] + REAL(1.0))/REAL(2.0);
-        summ += Gauss76Wt[i] * ellipsoid_kernel(q, b, a, zi);
+        //const real cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
+        const real cos_alpha = REAL(0.5)*(Gauss76Z[i] + REAL(1.0));
+        total += Gauss76Wt[i] * _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha);
     }
-    //const real form = (uplim-lolim)/2.0*summ;
-    const real form = REAL(0.5)*summ
-    const real s = (sld - sld_solvent) * form_volume(a, b);
+    //const real form = (upper-lower)/2*total;
+    const real form = REAL(0.5)*total;
+    const real s = (sld - solvent_sld) * form_volume(rpolar, requatorial);
     return REAL(1.0e-4) * form * s * s;
 }
@@ -68,6 +44,6 @@
     real sld,
     real solvent_sld,
-    real a,
-    real b,
+    real rpolar,
+    real requatorial,
     real theta,
     real phi)
@@ -77,7 +53,7 @@
     const real q = sqrt(qx*qx + qy*qy);
     SINCOS(theta*M_PI_180, sn, cn);
-    const real cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
-    const real form = ellipsoid_kernel(q, b, a, cos_val);
-    const real s = (sld - solvent_sld) * form_volume(a, b);
+    const real cos_alpha = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
+    const real form = _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha);
+    const real s = (sld - solvent_sld) * form_volume(rpolar, requatorial);
 
     return REAL(1.0e-4) * form * s * s;
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/ellipsoid.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,170 @@
+# ellipsoid model
+# Note: model title and parameter table are inserted automatically
+r"""
+The form factor is normalized by the particle volume.
+
+Definition
+----------
+
+The output of the 2D scattering intensity function for oriented ellipsoids
+is given by (Feigin, 1987)
+
+.. math::
+
+    P(q,\alpha) = \frac{\text{scale}}{V} f^2(q) + \text{background}
+
+where
+
+.. math::
+
+    f(q) = \frac{3 (\Delta rho)) V (\sin[qr(R_p,R_e,\alpha)] \
+                 - \cos[qr(R_p,R_e,\alpha)])}{[qr(R_q,R_b,\alpha)]^3}
+
+and
+
+.. math::
+
+    r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2}
+
+
+$\alpha$ is the angle between the axis of the ellipsoid and $\vec q$,
+$V$ is the volume of the ellipsoid, $R_p$ is the polar radius along the
+rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular
+to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the
+scattering length density difference between the scatterer and the solvent.
+
+To provide easy access to the orientation of the ellipsoid, we define
+the rotation axis of the ellipsoid using two angles $\theta$ and $\phi$.
+These angles are defined in Figure :num:`figure #cylinder-orientation`.
+For the ellipsoid, $\theta$ is the angle between the rotational axis
+and the $z$-axis.
+
+NB: The 2nd virial coefficient of the solid ellipsoid is calculated based
+on the $R_p$ and $R_e$ values, and used as the effective radius for
+$S(Q)$ when $P(Q) \dot S(Q)$ is applied.
+
+.. figure:: img/image121.JPG
+
+    The output of the 1D scattering intensity function for randomly oriented
+    ellipsoids given by the equation above.
+
+
+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.
+
+.. figure:: img/image122.JPG
+
+    The angles for oriented ellipsoid.
+
+Validation
+----------
+
+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 ellipsoid-comparison-1d` below shows a comparison of
+the 1D output of our model and the output of the NIST software.
+
+.. _ellipsoid-comparison-1d:
+
+.. figure:: img/image123.JPG
+
+    Comparison of the SasView scattering intensity for an ellipsoid
+    with the output of the NIST SANS analysis software.  The parameters
+    were set to: *scale=1.0*, *a=20 |Ang|*, *b=400 |Ang|*,
+    *sld-solvent_sld=3e-6 |Ang^-2|*, and *background=0.01 |cm^-1|*.
+
+Averaging over a distribution of orientation is done by evaluating the
+equation above. Since we have no other software to compare the
+implementation of the intensity for fully oriented ellipsoids, we can
+compare the result of averaging our 2D output using a uniform distribution
+$p(\theta,\phi) = 1.0$.  Figure :num:`figure #ellipsoid-comparison-2d`
+shows the result of such a cross-check.
+
+.. _ellipsoid-comparison-2d:
+
+.. figure:: img/image124.jpg
+
+    Comparison of the intensity for uniformly distributed ellipsoids
+    calculated from our 2D model and the intensity from the NIST SANS
+    analysis software. The parameters used were: *scale=1.0*,
+    *a=20 |Ang|*, *b=400 |Ang|*, *sld-solvent_sld=3e-6 |Ang^-2|*,
+    and *background=0.0 |cm^-1|*.
+
+The discrepancy above *q=0.3 |cm^-1|* is due to the way the form factors
+are calculated in the c-library provided by NIST. A numerical integration
+has to be performed to obtain *P(q)* for randomly oriented particles.
+The NIST software performs that integration with a 76-point Gaussian
+quadrature rule, which will become imprecise at high q where the amplitude
+varies quickly as a function of $q$. The SasView result shown has been
+obtained by summing over 501 equidistant points. Our result was found
+to be stable over the range of *q* shown for a number of points higher
+than 500.
+
+REFERENCE
+
+L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum,
+New York, 1987.
+"""
+
+from numpy import pi, inf
+
+name = "ellipsoid"
+title = "Ellipsoid of revolution with uniform scattering length density."
+
+description = """\
+P(q.alpha)= scale*f(q)^2 + background, where f(q)= 3*(sld
+		- solvent_sld)*V*[sin(q*r(Rp,Re,alpha))
+		-q*r*cos(qr(Rp,Re,alpha))]
+		/[qr(Rp,Re,alpha)]^3"
+
+     r(Rp,Re,alpha)= [Re^(2)*(sin(alpha))^2
+		+ Rp^(2)*(cos(alpha))^2]^(1/2)
+
+		sld: SLD of the ellipsoid
+		solvent_sld: SLD of the solvent
+		V: volume of the ellipsoid
+		Rp: polar radius of the ellipsoid
+		Re: equatorial radius of the ellipsoid
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "",
+      "Ellipsoid scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "",
+      "Solvent scattering length density" ],
+    [ "rpolar", "Ang",  20, [0, inf], "volume",
+      "Polar radius" ],
+    [ "requatorial", "Ang",  400, [0, inf], "volume",
+      "Equatorial radius" ],
+    [ "theta", "degrees", 60, [-inf, inf], "orientation",
+      "In plane angle" ],
+    [ "phi", "degrees", 60, [-inf, inf], "orientation",
+      "Out of plane angle" ],
+    ]
+
+source = [ "lib/J1.c", "lib/gauss76.c", "ellipsoid.c"]
+
+def ER(rpolar, requatorial):
+    import numpy as np
+
+    ee = np.empty_like(rpolar)
+    idx = rpolar > requatorial
+    ee[idx] = (rpolar[idx]**2 - requatorial[idx]**2)/rpolar[idx]**2
+    idx = rpolar < requatorial
+    ee[idx] = (requatorial[idx]**2 - rpolar[idx]**2)/requatorial[idx]**2
+    idx = rpolar == requatorial
+    ee[idx] = 2*rpolar[idx]
+    valid = (rpolar*requatorial != 0)
+    bd = 1.0-ee[valid]
+    e1 = np.sqrt(ee[valid])
+    b1 = 1.0 + np.arcsin(e1)/(e1*np.sqrt(bd))
+    bL = (1.0+e1)/(1.0-e1)
+    b2 = 1.0 + bd/2/e1*np.log(bL)
+    delta = 0.75*b1*b2
+
+    ddd = np.zeros_like(rpolar)
+    ddd[valid] = 2.0*(delta+1.0)*rpolar*requatorial**2
+    return 0.5*ddd**(1.0/3.0)
Index: sasmodels/models/lamellar.py
===================================================================
--- sasmodels/models/lamellar.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/lamellar.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,95 @@
+r"""
+Polydispersity in the bilayer thickness can be applied from the GUI.
+
+Definition
+----------
+
+The scattering intensity *I(q)* is
+
+.. math::
+
+    I(q) = 2\pi\frac{P(q){\delta q^2}
+
+
+The form factor is
+
+.. math::
+
+    P(q) = \frac{2\Delta\rho^2}{q^2}(1-cos(q\delta)
+
+
+where |delta| = bilayer thickness.
+
+The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 q_y^2}
+
+
+
+Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
+(Kline, 2006).
+
+.. figure:: img/image135.jpg
+
+    1D plot using the default values (w/1000 data point).
+
+
+Reference
+---------
+
+F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
+
+also in J. Phys. Chem. B, 105, (2001) 11081-11088
+
+
+"""
+
+from numpy import pi, inf
+
+name = "lamellar"
+title = "Lyotropic lamellar phase with uniform SLD and random distribution"
+description = """\
+	[Dilute Lamellar Form Factor](from a lyotropic lamellar phase)
+		I(q)= 2*pi*P(q)/(delta *q^(2)), where
+		P(q)=2*(contrast/q)^(2)*(1-cos(q*delta))^(2))
+		thickness = layer thickness
+		sld = layer scattering length density
+		sld_solvent = solvent scattering length density
+		background = incoherent background
+		scale = scale factor
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "sld", "1e-6/Ang^2", 1, [-inf,inf], "",
+      "Layer scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 6, [-inf,inf], "",
+      "Solvent scattering length density" ],
+    [ "thickness", "Ang",  50, [0, inf], "volume",
+      "Bilayer thickness" ],
+    ]
+
+
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return REAL(1.0);
+    """
+
+Iq = """
+    const real sub = sld - solvent_sld;
+    const real qsq = q*q;
+    return REAL(4.0e-4)*M_PI*sub*sub/qsq * (REAL(1.0)-cos(q*thickness))
+        / (thickness*qsq);
+    """
+
+Iqxy = """
+    // never called since no orientation or magnetic parameters.
+    return REAL(-1.0);
+    """
+
+# ER defaults to 0.0
+# VR defaults to 1.0
Index: smodels/models/lib/cylkernel.c
===================================================================
--- sasmodels/models/lib/cylkernel.c	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ 	(revision )
@@ -1,19 +1,0 @@
-// Compute the form factor for a cylinder
-//     F^2(q) sin a
-// given
-//     F(q) = sin(Q L/2 cos a)/(Q L/2 cos a) * 2 J1(Q R sin a) / (Q R sin a)
-// q is the q-value for the calculation (1/A)
-// radius is the radius of the cylinder (A)
-// halflength is the HALF-LENGTH of the cylinder = L/2 (A)
-real CylKernel(real q, real radius, real halflength, real theta);
-real CylKernel(real q, real radius, real halflength, real theta)
-{
-    real sn, cn;
-    SINCOS(theta, sn, cn);
-    const real besarg = q*radius*sn;
-    const real siarg = q*halflength*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)*sn*bj*bj*si*si;
-}
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/sphere.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,49 @@
+r"""
+
+"""
+
+from numpy import pi, inf
+
+name = "sphere"
+title = "Spheres with uniform SLD"
+description = """\
+	[Sphere Form Factor]
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "sld", "1e-6/Ang^2", 1, [-inf,inf], "",
+      "Layer scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 6, [-inf,inf], "",
+      "Solvent scattering length density" ],
+    [ "radius", "Ang",  50, [0, inf], "volume",
+      "Sphere radius" ],
+    ]
+
+
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return REAL(1.333333333333333)*M_PI*radius*radius*radius;
+    """
+
+Iq = """
+    real sn, cn;
+    const real qr = q*radius;
+    SINCOS(qr, sn, cn);
+    const real bes = (qr==REAL(0.0) ? REAL(1.0) : (REAL(3.0)*(sn-qr*cn))/(qr*qr*qr));
+    const real f = bes * (sld - solvent_sld) * form_volume(radius);
+    return REAL(1.0e-4)*f*f;
+    """
+
+
+Iqxy = """
+    // never called since no orientation or magnetic parameters.
+    return REAL(-1.0);
+    """
+
+def ER(radius):
+    return radius
+
+# VR defaults to 1.0
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,85 @@
+real form_volume(real req_minor, real req_major, real rpolar);
+real Iq(real q, real sld, real solvent_sld,
+    real req_minor, real req_major, real rpolar);
+real Iqxy(real qx, real qy, real sld, real solvent_sld,
+    real req_minor, real req_major, real rpolar, real theta, real phi, real psi);
+
+real form_volume(real req_minor, real req_major, real rpolar)
+{
+    return REAL(1.333333333333333)*M_PI*req_minor*req_major*rpolar;
+}
+
+real Iq(real q,
+    real sld,
+    real solvent_sld,
+    real req_minor,
+    real req_major,
+    real rpolar)
+{
+    // if (req_minor > req_major || req_major > rpolar) return REAL(-1.0);  // Exclude invalid region
+
+    real sn, cn;
+    real st, ct;
+    //const real lower = REAL(0.0);
+    //const real upper = REAL(1.0);
+    real outer = REAL(0.0);
+    for (int i=0;i<76;i++) {
+        //const real cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
+        const real x = REAL(0.5)*(Gauss76Z[i] + REAL(1.0));
+        SINCOS(M_PI_2*x, sn, cn);
+        const real acosx2 = req_minor*req_minor*cn*cn;
+        const real bsinx2 = req_major*req_major*sn*sn;
+        const real c2 = rpolar*rpolar;
+
+        real inner = REAL(0.0);
+        for (int j=0;j<76;j++) {
+            const real y = REAL(0.5)*(Gauss76Z[j] + REAL(1.0));
+            const real t = q*sqrt(acosx2 + bsinx2*(REAL(1.0)-y*y) + c2*y*y);
+            SINCOS(t, st, ct);
+            const real fq = ( t==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(st-t*ct)/(t*t*t) );
+            inner += Gauss76Wt[j] * fq * fq ;
+        }
+        outer += Gauss76Wt[i] * REAL(0.5) * inner;
+    }
+    //const real fq2 = (upper-lower)/2*outer;
+    const real fq2 = REAL(0.5)*outer;
+    const real s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);
+    return REAL(1.0e-4) * fq2 * s * s;
+}
+
+real Iqxy(real qx, real qy,
+    real sld,
+    real solvent_sld,
+    real req_minor,
+    real req_major,
+    real rpolar,
+    real theta,
+    real phi,
+    real psi)
+{
+    // if (req_minor > req_major || req_major > rpolar) return REAL(-1.0);  // Exclude invalid region
+
+    real stheta, ctheta;
+    real sphi, cphi;
+    real spsi, cpsi;
+    real st, ct;
+
+    const real q = sqrt(qx*qx + qy*qy);
+    const real qxhat = qx/q;
+    const real qyhat = qy/q;
+    SINCOS(theta*M_PI_180, stheta, ctheta);
+    SINCOS(phi*M_PI_180, sphi, cphi);
+    SINCOS(psi*M_PI_180, spsi, cpsi);
+    const real calpha = ctheta*cphi*qxhat + stheta*qyhat;
+    const real cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat;
+    const real cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat;
+    const real t = q*sqrt(req_minor*req_minor*cnu*cnu
+                          + req_major*req_major*cmu*cmu
+                          + rpolar*rpolar*calpha*calpha);
+    SINCOS(t, st, ct);
+    const real fq = ( t==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(st-t*ct)/(t*t*t) );
+    const real s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);
+
+    return REAL(1.0e-4) * fq * fq * s * s;
+}
+
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -0,0 +1,130 @@
+# triaxial ellipsoid model
+# Note: model title and parameter table are inserted automatically
+r"""
+All three axes are of different lengths with $R_a \le R_b <= R_c$
+**Users should maintain this inequality for all calculations**.
+
+.. math::
+
+    P(q) = \text{scale} V \left< f^2(q) \right> + \text{background}
+
+where the volume $V = 4/3 \pi R_a R_b R_c$, and the averaging
+$\left< \cdots \right>$ is applied over all orientations for 1D.
+
+.. figure:: img/image128.JPG
+
+    Ellipsoid schematic.
+
+The returned value is in units of |cm^-1|, on absolute scale.
+
+Definition
+----------
+
+The form factor calculated is
+
+.. math::
+
+    P(q) = \frac{\text{scale}}{V}\int_0^1\int_0^1 \
+        \Phi^2(qR_a^2\cos^2( \pi x/2) + qR_b^2\sin^2(\pi y/2)(1-y^2) + c^2y^2) \
+        dx dy
+
+where
+
+.. math::
+
+    \Phi(u) = 3 u^{-3} (\sin u - u \cos u)
+
+To provide easy access to the orientation of the triaxial ellipsoid,
+we define the axis of the cylinder using the angles $\theta$, $\phi$
+and $\psi$. These angles are defined on
+Figure :num:`figure #triaxial-ellipsoid-angles`.
+The angle $\psi$ is the rotational angle around its own $c$ axis
+against the $q$ plane. For example, $\psi = 0$ when the
+$a$ axis is parallel to the $x$ axis of the detector.
+
+.. _triaxial-ellipsoid-angles:
+
+.. figure:: img/image132.JPG
+
+    The angles for oriented ellipsoid.
+
+The radius-of-gyration for this system is  $Rg^2 = (R_a R_b R_c)^2/5$.
+
+The contrast is defined as SLD(ellipsoid) - SLD(solvent).  In the
+parameters, *a* is the minor equatorial radius, *b* is the major
+equatorial radius, and c is the polar radius of the ellipsoid.
+
+NB: The 2nd virial coefficient of the triaxial solid ellipsoid is
+calculated based on the polar radius $R_p = R_c$ and equatorial
+radius $R_e = \sqrt{R_a R_b}$, and used as the effective radius for
+$S(Q)$ when $P(Q) \dot S(Q)$ is applied.
+
+.. figure:: img/image130.JPG
+
+    1D plot using the default values (w/1000 data point).
+
+Validation
+----------
+
+Validation of our code was done by comparing the output of the
+1D calculation to the angular average of the output of 2D calculation
+over all possible angles.
+Figure :num:`figure #triaxial-ellipsoid-compare` shows the comparison where
+the solid dot refers to averaged 2D while the line represents the
+result of 1D calculation (for 2D averaging, 76, 180, and 76 points
+are taken for the angles of $\theta$, $\phi$, and $\psi$ respectively).
+
+.. _triaxial-ellipsoid-compare:
+
+.. figure:: img/image131.GIF
+
+    Comparison between 1D and averaged 2D.
+
+Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
+(Kline, 2006)
+
+REFERENCE
+
+L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum,
+New York, 1987.
+"""
+
+from numpy import pi, inf
+
+name = "triaxial_ellipsoid"
+title = "Ellipsoid of uniform scattering length density with three independent axes."
+
+description = """\
+Note: During fitting ensure that the inequality ra<rb<rc is not
+	violated. Otherwise the calculation will
+	not be correct.
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "",
+      "Ellipsoid scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "",
+      "Solvent scattering length density" ],
+    [ "req_minor", "Ang",  20, [0, inf], "volume",
+      "Minor equitorial radius" ],
+    [ "req_major", "Ang",  400, [0, inf], "volume",
+      "Major equatorial radius" ],
+    [ "rpolar", "Ang",  10, [0, inf], "volume",
+      "Polar radius" ],
+    [ "theta", "degrees", 60, [-inf, inf], "orientation",
+      "In plane angle" ],
+    [ "phi", "degrees", 60, [-inf, inf], "orientation",
+      "Out of plane angle" ],
+    [ "psi", "degrees", 60, [-inf, inf], "orientation",
+      "Out of plane angle" ],
+    ]
+
+source = [ "lib/J1.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
+
+def ER(req_minor, req_major, rpolar):
+    import numpy as np
+    from .ellipsoid import ER as ellipsoid_ER
+    return ellipsoid_ER(rpolar, np.sqrt(req_minor*req_major))
+
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/sasview_model.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -6,9 +6,9 @@
 
 try:
-    import pyopencl
-    from .gen import opencl_model as load_model
-except ImportError:
+    from .gpu import load_model
+except ImportError,exc:
+    warnings.warn(str(exc))
     warnings.warn("OpenCL not available --- using ctypes instead")
-    from .gen import dll_model as load_model
+    from .dll import load_model
 
 
@@ -246,5 +246,9 @@
             # Check whether we have a list of ndarrays [qx,qy]
             qx, qy = qdist
-            return self.calculate_Iq(qx, qy)
+            partype = self._model.info['partype']
+            if not partype['orientation'] and not partype['magnetic']:
+                return self.calculate_Iq(np.sqrt(qx**2+qy**2))
+            else:
+                return self.calculate_Iq(qx, qy)
 
         elif isinstance(qdist, np.ndarray):
Index: sasmodels/weights.py
===================================================================
--- sasmodels/weights.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/weights.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
@@ -42,6 +42,9 @@
         """
         sigma = self.width * center if relative else self.width
-        if sigma == 0:
-            return np.array([center], 'd'), np.array([1.], 'd')
+        if sigma == 0 or self.npts < 2:
+            if lb <= center <= ub:
+                return np.array([center], 'd'), np.array([1.], 'd')
+            else:
+                return np.array([], 'd'), np.array([], 'd')
         return self._weights(center, sigma, lb, ub)
 
