Index: Kernel/Kernel-Cylinder_f.cpp
===================================================================
--- Kernel/Kernel-Cylinder_f.cpp	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ Kernel/Kernel-Cylinder_f.cpp	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -0,0 +1,113 @@
+#ifdef __OPENCL_VERSION__
+# define USE_OPENCL
+#endif
+
+// If opencl is not available
+#ifndef USE_OPENCL
+#  include <math.h>
+#  define REAL(x) (x)
+#  define real double
+#  define global
+#  define local
+#  define kernel extern "C"
+#  include "NR_BessJ1.cpp"
+#endif
+
+#ifndef M_PI
+#  define M_PI REAL(3.141592653589793)
+#endif
+#define RADIANS REAL(0.017453292519943295)
+#define PI_E8 REAL(3.141592653589793e8)
+real f(real qx, real qy, real sub, real rr, real h, real scale, real cyl_theta, real cyl_phi);
+real f(real qx, real qy, real sub, real rr, real h, real scale, real cyl_theta, real cyl_phi)
+{
+    const real qq = sqrt(qx*qx+qy*qy);
+
+    const real theta = cyl_theta*RADIANS;
+    const real cos_val = cos(theta)*cos(cyl_phi*RADIANS)*(qx/qq) + sin(theta)*(qy/qq);
+
+    real alpha = acos(cos_val);
+    if(alpha == REAL(0.0)) {
+        alpha = REAL(1.0e-26);
+    }
+    const real sin_alpha = sin(alpha);
+    const real besarg = qq*rr*sin_alpha;
+    const real siarg = qq*h/2*cos(alpha);
+
+    real be, si;
+    if (besarg == REAL(0.0)) {
+        be = sin_alpha;
+    } else {
+        const real bj = NR_BessJ1(besarg)/besarg;
+        be = REAL(4.0)*bj*bj*sin_alpha;
+    }
+    if (siarg == REAL(0.0)) {
+        si = REAL(1.0);
+    } else {
+        si = sin(siarg)/siarg;
+    }
+    const real form = be*si*si/sin_alpha;
+    return PI_E8*scale*form*sub*sub*rr*rr*rr*rr*h*h;
+}
+
+kernel void CylinderKernel(global const real *qx, global const real *qy, global real *result,
+#ifdef USE_OPENCL
+    global real *loops_g,
+#else
+    const int Nq,
+#endif
+    local real *loops, const real cutoff,
+    const real scale, const real background,
+    const real sub,
+    const int Nradius, const int Nlength, const int Ntheta, const int Nphi)
+{
+#ifdef USE_OPENCL
+    // copy loops info to local memory
+    event_t e = async_work_group_copy(loops, loops_g, (Nradius+Nlength+Ntheta+Nphi)*2, 0);
+    wait_group_events(1, &e);
+
+    int i = get_global_id(0);
+    int count = get_global_size(0);
+
+    if(i < count)
+#else
+    #pragma omp parallel for
+    for (int i=0; i<Nq; i++)
+#endif
+    {
+        const real qxi=qx[i];
+        const real qyi=qy[i];
+        real ret=REAL(0.0), norm=REAL(0.0), norm_vol=REAL(0.0), vol=REAL(0.0);
+        for (int ri=0; ri < Nradius; ri++) {
+            const real rv = loops[2*ri];
+            const real rw = loops[2*ri+1];
+            for (int li=0; li < Nlength; li++) {
+                const real lv = loops[2*(li+Nradius)];
+                const real lw = loops[2*(li+Nradius)+1];
+                vol += rw*lw*rv*rv*lv*Ntheta*Nphi;
+               norm_vol += rw*lw*Ntheta*Nphi;
+                for (int thi=0; thi < Ntheta; thi++) {
+                    const real thv = loops[2*(thi+Nradius+Nlength)];
+                    const real thw = loops[2*(thi+Nradius+Nlength)+1];
+                    // #pragma unroll
+                    for (int phi=0; phi < Nphi; phi++) {
+                        const real phv = loops[2*(phi+Nradius+Nlength+Ntheta)];
+                        const real phw = loops[2*(phi+Nradius+Nlength+Ntheta)+1];
+
+                        const real weight = rw*lw*thw*phw;
+                        //ret += qxi + qyi + sub + rv + lv + weight + thv + phv;
+                        if (weight > cutoff) {
+                            ret += f(qxi, qyi, sub, rv, lv, weight, thv, phv);
+                            norm += weight;
+                        }
+                    }
+                }
+            }
+        }
+        if (Ntheta>1) norm = norm/(M_PI/2);
+        if (vol != REAL(0.0) && norm_vol != REAL(0.0)) {
+            ret *= norm_vol/vol;
+        }
+        result[i] = scale*ret/norm+background;
+    }
+}
Index: Models/code_capcyl.py
===================================================================
--- Models/code_capcyl.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ Models/code_capcyl.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -12,8 +12,8 @@
 
 class GpuCapCylinder(object):
-    PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'length':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0,
+    PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'len_cyl':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0,
              'theta':0, 'phi':0}
 
-    PD_PARS = ['rad_cyl', 'length', 'rad_cap', 'theta', 'phi']
+    PD_PARS = ['rad_cyl', 'len_cyl', 'rad_cap', 'theta', 'phi']
 
     def __init__(self, qx, qy, dtype='float32'):
@@ -39,5 +39,6 @@
         self.res[:] = 0
         cl.enqueue_copy(queue, self.res_b, self.res)
-        rad_cyl,length,rad_cap,theta,phi = \
+
+        rad_cyl,len_cyl,rad_cap,theta,phi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
              for base in GpuCapCylinder.PD_PARS]
@@ -45,5 +46,5 @@
         rad_cyl.value, rad_cyl.weight = rad_cyl.get_weights(pars['rad_cyl'], 0, 10000, True)
         rad_cap.value, rad_cap.weight = rad_cap.get_weights(pars['rad_cap'], 0, 10000, True)
-        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
+        len_cyl.value, len_cyl.weight = len_cyl.get_weights(pars['len_cyl'], 0, 10000, True)
         theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False)
         phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False)
@@ -56,19 +57,24 @@
         for i in xrange(len(rad_cyl.weight)):
             for m in xrange(len(rad_cap.weight)):
-                for j in xrange(len(length.weight)):
+                for j in xrange(len(len_cyl.weight)):
+
+                    hDist = -1.0*sqrt(fabs(rad_cap.value[m]*rad_cap.value[m]-rad_cyl.value[i]*rad_cyl.value[i]))
+                    vol_i = 4.0*atan(1.0)*rad_cyl.value[i]*rad_cyl.value[i]*len_cyl.value[j]+2.0*4.0*atan(1.0)/3.0\
+                                    *((rad_cap.value[m]-hDist)*(rad_cap.value[m]-hDist)*(2*rad_cap.value[m]+hDist))
+                    vol += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]*vol_i
+                    norm_vol += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]
+
                     for k in xrange(len(theta.weight)):
                         for l in xrange(len(phi.weight)):
-                            hDist = -1.0*sqrt(fabs(rad_cap.value[m]*rad_cap.value[m]-rad_cyl.value[i]*rad_cyl.value[i]))
-                            vol_i = 4.0*atan(1.0)*rad_cyl.value[i]*rad_cyl.value[i]*length.value[j]+2.0*4.0*atan(1.0)/3.0\
-                                            *((rad_cap.value[m]-hDist)*(rad_cap.value[m]-hDist)*(2*rad_cap.value[m]+hDist))
+
                             self.prg.CapCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                        real(vol_i), real(hDist), real(rad_cyl.value[i]), real(rad_cap.value[m]), real(length.value[j]),
+                                        real(vol_i), real(hDist), real(rad_cyl.value[i]), real(rad_cap.value[m]), real(len_cyl.value[j]),
                                         real(theta.value[k]), real(phi.value[l]), real(sub), real(pars['scale']),
                                         real(phi.weight[l]), real(theta.weight[k]), real(rad_cap.weight[m]),
-                                        real(rad_cyl.weight[i]), real(length.weight[j]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size))
+                                        real(rad_cyl.weight[i]), real(len_cyl.weight[j]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size))
 
-                            vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*vol_i
-                            norm_vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]
-                            norm += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*theta.weight[k]*phi.weight[l]
+                            norm += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]*theta.weight[k]*phi.weight[l]
+
+
 
         #if size > 1:
Index: Models/code_coreshellcyl.py
===================================================================
--- Models/code_coreshellcyl.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ Models/code_coreshellcyl.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -48,7 +48,12 @@
         for r in xrange(len(radius.weight)):
             for l in xrange(len(length.weight)):
-                for at in xrange(len(axis_theta.weight)):
-                    for p in xrange(len(axis_phi.weight)):
-                        for th in xrange(len(thickness.weight)):
+                for th in xrange(len(thickness.weight)):
+
+                    vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\
+                                   *(length.value[l]+2.0*thickness.value[th])
+                    norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th]
+
+                    for at in xrange(len(axis_theta.weight)):
+                        for p in xrange(len(axis_phi.weight)):
 
                             self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
@@ -59,7 +64,5 @@
                                     real(pars['shell_sld']), real(pars['solvent_sld']),np.uint32(size),
                                     np.uint32(self.qx.size))
-                            vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\
-                                   *(length.value[l]+2.0*thickness.value[th])
-                            norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th]
+
                             norm += radius.weight[r]*length.weight[l]*thickness.weight[th]*axis_theta.weight[at]\
                                     *axis_phi.weight[p]
Index: Models/code_cylinder.py
===================================================================
--- Models/code_cylinder.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ Models/code_cylinder.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -6,5 +6,5 @@
 
 from weights import GaussianDispersion
-from sasmodel import card, set_precision, set_precision_1d
+from sasmodel import card, set_precision, set_precision_1d, tic, toc
 
 
@@ -33,4 +33,5 @@
     def eval(self, pars):
 
+        tic()
         _ctx,queue = card()
         self.res[:] = 0
@@ -55,4 +56,8 @@
         for i in xrange(len(radius.weight)):
             for j in xrange(len(length.weight)):
+
+                vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j]
+                norm_vol += radius.weight[i]*length.weight[j]
+
                 for k in xrange(len(cyl_theta.weight)):
                     for l in xrange(len(cyl_phi.weight)):
@@ -63,7 +68,6 @@
                                            np.uint32(self.qx.size), np.uint32(size))
 
-                        vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j]
-                        norm_vol += radius.weight[i]*length.weight[j]
                         norm += radius.weight[i]*length.weight[j]*cyl_theta.weight[k]*cyl_phi.weight[l]
+
 
        # if size > 1:
@@ -74,4 +78,5 @@
             sum *= norm_vol/vol
 
+        print toc()*1000, self.qx.shape[0]
         return sum/norm+pars['background']
 
Index: Models/code_cylinder_f.py
===================================================================
--- Models/code_cylinder_f.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ Models/code_cylinder_f.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -0,0 +1,147 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+import ctypes
+from ctypes import c_int, c_double, c_void_p
+import numpy as np
+import pyopencl as cl
+from pyopencl import mem_flags as mf
+
+from weights import GaussianDispersion
+from sasmodel import card, set_precision, set_precision_1d, tic, toc
+
+
+class GpuCylinder(object):
+    PARS = {
+        'scale':1,'radius':1,'length':1,'sldCyl':1e-6,'sldSolv':0,'background':0,
+        'cyl_theta':0,'cyl_phi':0,
+    }
+    PD_PARS = ['radius', 'length', 'cyl_theta', 'cyl_phi']
+
+    def __init__(self, qx, qy, dtype='float32', cutoff=1e-5):
+
+        #create context, queue, and build program
+        ctx,_queue = card()
+        bessel = open('Kernel/NR_BessJ1.cpp').read()
+        kernel = open('Kernel/Kernel-Cylinder_f.cpp').read()
+        src, qx, qy = set_precision("\n".join((bessel,kernel)), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+        self.cutoff = cutoff
+
+        #buffers
+        self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
+        self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
+        self.res_b = cl.Buffer(ctx, cl.mem_flags.READ_WRITE, self.qx.nbytes)
+        self.res = np.empty_like(self.qx)
+
+    def eval(self, pars):
+        tic()
+
+        ctx,queue = card()
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
+        loops, loop_lengths = make_loops(pars, dtype=self.qx.dtype)
+        loops_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops)
+        loops_l = cl.LocalMemory(len(loops.data))
+        self.prg.CylinderKernel(queue, self.qx.shape, None,
+            self.qx_b, self.qy_b, self.res_b,
+            loops_b, loops_l, real(self.cutoff),
+            real(pars['scale']), real(pars['background']),
+            real(pars['sldCyl']-pars['sldSolv']),
+            *[np.uint32(pn) for pn in loop_lengths])
+        cl.enqueue_copy(queue, self.res, self.res_b)
+        print toc()*1000, self.qx.shape[0]
+
+        return self.res
+
+class CpuCylinder(GpuCylinder): #inherit parameters only
+    def __init__(self, qx, qy, dtype='float32', cutoff=1e-5):
+        self.qx, self.qy = [np.ascontiguousarray(v,'d') for v in qx,qy]
+        self.cutoff = cutoff
+        self.res = np.empty_like(self.qx)
+        self.dll = ctypes.CDLL('Kernel/cylinder.so')
+        self.fn = self.dll['CylinderKernel']
+        self.fn.argtypes = [
+            c_void_p, c_void_p, c_void_p, c_int,
+            c_void_p, c_double, c_double, c_double, c_double,
+            c_int, c_int, c_int, c_int
+        ]
+    def eval(self, pars):
+        loops, loop_lengths = make_loops(pars, dtype=self.qx.dtype)
+        self.fn(self.qx.ctypes.data, self.qy.ctypes.data, self.res.ctypes.data, len(self.qx),
+            loops.ctypes.data, self.cutoff, pars['scale'], pars['background'],
+            pars['sldCyl']-pars['sldSolv'], *loop_lengths)
+
+        return self.res
+
+def make_loops(pars, dtype='double'):
+    # 0.2 ms on sparkle to form the final loops
+    radius, length, theta, phi = \
+       [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+         for base in GpuCylinder.PD_PARS]
+    parts = [
+        radius.get_weights(pars['radius'], 0, 10000, True),
+        length.get_weights(pars['length'], 0, 10000, True),
+        theta.get_weights(pars['cyl_theta'], -np.inf, np.inf, False),
+        theta.get_weights(pars['cyl_phi'], -np.inf, np.inf, False),
+        ]
+    # Make sure that weights are normalized to peaks at 1 so that
+    # the tolerance term can be used properly on truncated distributions
+    loops = np.hstack((v,w/w.max()) for v,w in parts)
+    #loops = np.hstack(parts)
+    loops = np.ascontiguousarray(loops.T, dtype).flatten()
+    return loops, [len(p[0]) for p in parts]
+
+class OneDGpuCylinder(object):
+    PARS = {
+        'scale':1,'radius':1,'length':1,'sldCyl':1e-6,'sldSolv':0,'background':0,
+        'bolim':0, 'uplim':90
+    }
+    PD_PARS = ['radius', 'length']
+
+    def __init__(self, q, dtype='float32'):
+
+        #create context, queue, and build program
+        ctx,_queue = card()
+        trala = open('Kernel/NR_BessJ1.cpp').read()+"\n"+open('Kernel/OneDCyl_Kfun.cpp').read()+"\n"+open('Kernel/Kernel-OneDCylinder.cpp').read()
+        src, self.q = set_precision_1d(trala, q, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+
+        #buffers
+        mf = cl.mem_flags
+        self.q_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.q)
+        self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, q.nbytes)
+        self.res = np.empty_like(self.q)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        radius, length = \
+            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+             for base in OneDGpuCylinder.PD_PARS]
+
+        #Get the weights for each
+        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True)
+        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
+
+        #Perform the computation, with all weight points
+        sum, norm, vol = 0.0, 0.0, 0.0,
+        sub = pars['sldCyl'] - pars['sldSolv']
+
+        real = np.float32 if self.q.dtype == np.dtype('float32') else np.float64
+        #Loop over radius, length, theta, phi weight points
+        for r in xrange(len(radius.weight)):
+            for l in xrange(len(length.weight)):
+                        self.prg.OneDCylKernel(queue, self.q.shape, None, self.q_b, self.res_b, real(sub),
+                                           real(length.value[l]), real(radius.value[r]), real(pars['scale']),
+                                           np.uint32(self.q.size), real(pars['uplim']), real(pars['bolim']))
+                        cl.enqueue_copy(queue, self.res, self.res_b)
+                        sum += radius.weight[r]*length.weight[l]*self.res*pow(radius.value[r],2)*length.value[l]
+                        vol += radius.weight[r]*length.weight[l] *pow(radius.value[r],2)*length.value[l]
+                        norm += radius.weight[r]*length.weight[l]
+
+        if vol != 0.0 and norm != 0.0:
+            sum *= norm/vol
+
+        return sum/norm + pars['background']
Index: Models/code_ellipse.py
===================================================================
--- Models/code_ellipse.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ Models/code_ellipse.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -53,4 +53,7 @@
             for j in xrange(len(radius_b.weight)):
                 #Average over theta distribution
+                vol += radius_a.weight[i]*radius_b.weight[j]*pow(radius_b.value[j], 2)*radius_a.value[i]
+                norm_vol += radius_a.weight[i]*radius_b.weight[j]
+
                 for k in xrange(len(axis_theta.weight)):
                     #Average over phi distribution
@@ -64,7 +67,6 @@
                                         np.uint32(self.qx.size), np.uint32(len(axis_theta.weight)))
 
-                        vol += radius_a.weight[i]*radius_b.weight[j]*pow(radius_b.value[j], 2)*radius_a.value[i]
-                        norm_vol += radius_a.weight[i]*radius_b.weight[j]
                         norm += radius_a.weight[i]*radius_b.weight[j]*axis_theta.weight[k]*axis_phi.weight[l]
+
         # Averaging in theta needs an extra normalization
         # factor to account for the sin(theta) term in the
Index: Models/code_triaxialellipse.py
===================================================================
--- Models/code_triaxialellipse.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ Models/code_triaxialellipse.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -9,8 +9,8 @@
 
 class GpuTriEllipse:
-    PARS = {'scale':1, 'axisA':35, 'axisB':100, 'axisC':400, 'sldEll':1e-6, 'sldSolv':6.3e-6, 'background':0,
-            'theta':0, 'phi':0, 'psi':0}
+    PARS = {'scale':1, 'semi_axisA':35, 'semi_axisB':100, 'semi_axisC':400, 'sldEll':1e-6, 'sldSolv':6.3e-6, 'background':0,
+            'axis_theta':0, 'axis_phi':0, 'axis_psi':0}
 
-    PD_PARS = ['axisA', 'axisB', 'axisC', 'theta', 'phi', 'psi']
+    PD_PARS = ['semi_axisA', 'semi_axisB', 'semi_axisC', 'axis_theta', 'axis_phi', 'axis_psi']
 
     def __init__(self, qx, qy, dtype='float32'):
@@ -32,35 +32,38 @@
         self.res[:] = 0
         cl.enqueue_copy(queue, self.res_b, self.res)
-        axisA, axisB, axisC, theta, phi, psi = \
+        semi_axisA, semi_axisB, semi_axisC, axis_theta, axis_phi, axis_psi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
              for base in GpuTriEllipse.PD_PARS]
 
-        axisA.value, axisA.weight = axisA.get_weights(pars['axisA'], 0, 10000, True)
-        axisB.value, axisB.weight = axisB.get_weights(pars['axisB'], 0, 10000, True)
-        axisC.value, axisC.weight = axisC.get_weights(pars['axisC'], 0, 10000, True)
-        theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False)
-        phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False)
-        psi.value, psi.weight = psi.get_weights(pars['psi'], -90, 180, False)
+        semi_axisA.value, semi_axisA.weight = semi_axisA.get_weights(pars['semi_axisA'], 0, 10000, True)
+        semi_axisB.value, semi_axisB.weight = semi_axisB.get_weights(pars['semi_axisB'], 0, 10000, True)
+        semi_axisC.value, semi_axisC.weight = semi_axisC.get_weights(pars['semi_axisC'], 0, 10000, True)
+        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False)
+        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False)
+        axis_psi.value, axis_psi.weight = axis_psi.get_weights(pars['axis_psi'], -90, 180, False)
 
         sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
-        size = len(theta.weight)
+        size = len(axis_theta.weight)
         sub = pars['sldEll'] - pars['sldSolv']
 
         real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-        for a in xrange(len(axisA.weight)):
-            for b in xrange(len(axisB.weight)):
-                for c in xrange(len(axisC.weight)):
-                    for t in xrange(len(theta.weight)):
-                        for i in xrange(len(phi.weight)):
-                            for s in xrange(len(psi.weight)):
+        for a in xrange(len(semi_axisA.weight)):
+            for b in xrange(len(semi_axisB.weight)):
+                for c in xrange(len(semi_axisC.weight)):
+
+                    vol += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]*semi_axisA.value[a]*semi_axisB.value[b]*semi_axisC.value[c]
+                    norm_vol += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]
+
+                    for t in xrange(len(axis_theta.weight)):
+                        for i in xrange(len(axis_phi.weight)):
+                            for s in xrange(len(axis_psi.weight)):
                                 self.prg.TriaxialEllipseKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                            real(sub), real(pars['scale']), real(axisA.value[a]), real(axisB.value[b]),
-                                            real(axisC.value[c]), real(phi.value[i]), real(theta.value[t]), real(psi.value[s]),
-                                            real(axisA.weight[a]), real(axisB.weight[b]), real(axisC.weight[c]), real(psi.weight[s]),
-                                            real(phi.weight[i]), real(theta.weight[t]), np.uint32(self.qx.size), np.uint32(size))
+                                            real(sub), real(pars['scale']), real(semi_axisA.value[a]), real(semi_axisB.value[b]),
+                                            real(semi_axisC.value[c]), real(axis_phi.value[i]), real(axis_theta.value[t]), real(axis_psi.value[s]),
+                                            real(semi_axisA.weight[a]), real(semi_axisB.weight[b]), real(semi_axisC.weight[c]), real(axis_psi.weight[s]),
+                                            real(axis_phi.weight[i]), real(axis_theta.weight[t]), np.uint32(self.qx.size), np.uint32(size))
 
-                                vol += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*axisA.value[a]*axisB.value[b]*axisC.value[c]
-                                norm_vol += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]
-                                norm += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*theta.weight[t]*phi.weight[i]*psi.weight[s]
+                                norm += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]*axis_theta.weight[t]*axis_phi.weight[i]*axis_psi.weight[s]
+
 
       #  if size > 1:
Index: compare.py
===================================================================
--- compare.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ compare.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -2,17 +2,7 @@
 # -*- coding: utf-8 -*-
 
-import datetime
-
-from sasmodel import SasModel, load_data, set_beam_stop, plot_data
-
-
-TIC = None
-def tic():
-    global TIC
-    TIC = datetime.datetime.now()
-
-def toc():
-    now = datetime.datetime.now()
-    return (now-TIC).total_seconds()
+from sasmodel import SasModel, load_data, set_beam_stop, plot_data, tic, toc
+import numpy as np
+
 
 def sasview_model(modelname, **pars):
@@ -40,7 +30,35 @@
     return theory
 
-def cyl(N=1):
+def plot(data, cpumodel, gpumodel, N, pars):
+
+    model = SasModel(data, gpumodel, dtype='f', **pars)
+    tic()
+    for i in range(N):
+        #pars['scale'] = np.random.rand()
+        model.update()
+        gpu = model.theory()
+    gpu_time = toc()*1000./N
+    print "ocl t=%.1f ms"%gpu_time
+    tic()
+    for i in range(N):
+        cpu = sasview_eval(cpumodel, data)
+    cpu_time = toc()*1000./N
+    print "omp t=%.1f ms"%cpu_time
+
+    import matplotlib.pyplot as plt
+
+    print "gpu/cpu", max(gpu/cpu)
+    cpu *= max(gpu/cpu)
+    relerr = (gpu - cpu)/cpu
+    print "max(|(ocl-omp)/ocl|)", max(abs(relerr))
+
+
+    plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time)
+    plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time)
+    plt.subplot(133); plot_data(data, relerr); plt.title("relerr x 10^8"); plt.colorbar()
+    plt.show()
+
+def newlen(N):
     import sys
-    import matplotlib.pyplot as plt
 
     if len(sys.argv) > 1:
@@ -49,85 +67,42 @@
     set_beam_stop(data, 0.004)
 
+    return data, N
+
+def cyl(N=1):
+
     dtype = "float"
     pars = dict(
-        scale=.08, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1,
-        cyl_theta=0, cyl_phi=0,
+        scale=.003, radius=64.1, length=66.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1,
+        cyl_theta=90, cyl_phi=0,
         radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
         length_pd=0.1,length_pd_n=5, length_pd_nsigma=3,
         cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3,
-        cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3,
-        )
-
+        cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3)
+
+    from Models.code_cylinder_f import GpuCylinder as gpumodel
     model = sasview_model('Cylinder', **pars)
-    tic()
-    for i in range(N):
-        cpu = sasview_eval(model, data)
-    cpu_time = toc()*1000./N
-
-
-    from Models.code_cylinder import GpuCylinder
-    model = SasModel(data, GpuCylinder, dtype='f', **pars)
-    tic()
-    for i in range(N):
-        gpu = model.theory()
-    gpu_time = toc()*1000./N
-
-    print "gpu/cpu", max(gpu/cpu)
-    cpu *= max(gpu/cpu)
-    relerr = (gpu - cpu)/cpu
-    print "max(|(ocl-omp)/ocl|)", max(abs(relerr))
-    print "omp t=%.1f ms"%cpu_time
-    print "ocl t=%.1f ms"%gpu_time
-
-    plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time)
-    plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time)
-    plt.subplot(133); plot_data(data, relerr); plt.title("relerr x 10^8"); plt.colorbar()
-    plt.show()
-
-def ellipse(N=4):
-    import sys
-    import matplotlib.pyplot as plt
-
-    if len(sys.argv) > 1:
-        N = int(sys.argv[1])
-    data = load_data('DEC07106.DAT')
-    set_beam_stop(data, 0.004)
+    data, N = newlen(N)
+
+    plot(data, model, gpumodel, N, pars)
+
+
+def ellipse(N=1):
 
     pars = dict(scale=.027, radius_a=60, radius_b=180, sldEll=.297e-6, sldSolv=5.773e-6, background=4.9,
-                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,)
-
+                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)
-    tic()
-    for i in range(N):
-        cpu = sasview_eval(model, data)
-    cpu_time = toc()*1000./N
-
-    from Models.code_ellipse import GpuEllipse
-    model = SasModel(data, GpuEllipse, dtype='f', **pars)
-    tic()
-    for i in range(N):
-        gpu = model.theory()
-    gpu_time = toc()*1000./N
-
-    relerr = (gpu - cpu)/cpu
-    print "max(|(ocl-omp)/ocl|)", max(abs(relerr))
-    print "omp t=%.1f ms"%cpu_time
-    print "ocl t=%.1f ms"%gpu_time
-
-    plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time)
-    plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time)
-    plt.subplot(133); plot_data(data, 1e8*relerr); plt.title("relerr x 10^8"); plt.colorbar()
-    plt.show()
+
+    data, N = newlen(N)
+    plot(data, model, gpumodel, N, pars)
 
 def coreshell(N=1):
-    import sys
-    import matplotlib.pyplot as plt
-
-    if len(sys.argv) > 1:
-        N = int(sys.argv[1])
-    data = load_data('December/DEC07098.DAT')
-    set_beam_stop(data, 0.004)
+
+    data, N = newlen(N)
 
     pars = dict(scale= 1.77881e-06, radius=325, thickness=25, length=34.2709,
@@ -135,62 +110,98 @@
                  background=223.827, axis_theta=90, axis_phi=0,
                  axis_theta_pd=15.8,
-                 radius_pd=0.1, radius_pd_n=1, radius_pd_nsigma=0,
-                 length_pd=0.1, length_pd_n=1, length_pd_nsigma=0,
-                 thickness_pd=0.1, thickness_pd_n=1, thickness_pd_nsigma=0,
-                 axis_theta_pd_n=20, axis_theta_pd_nsigma=3,
-                 axis_phi_pd=0.0008748, axis_phi_pd_n=60, axis_phi_pd_nsigma=3,)
+                 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_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)
-    tic()
-    for i in range(N):
-        cpu = sasview_eval(model, data)
-    cpu_time = toc()*1000./N
-
-    from Models.code_coreshellcyl import GpuCoreShellCylinder
-    model = SasModel(data, GpuCoreShellCylinder, dtype='f', **pars)
-    tic()
-    for i in range(N):
-        gpu = model.theory()
-    gpu_time = toc()*1000./N
-
-    relerr = (gpu - cpu)/cpu
-    print "max(|(ocl-omp)/ocl|)", max(abs(relerr))
-    print "omp t=%.1f ms"%cpu_time
-    print "ocl t=%.1f ms"%gpu_time
-
-    plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time)
-    plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time)
-    plt.subplot(133); plot_data(data, 1e8*relerr); plt.title("relerr x 10^8"); plt.colorbar()
-    plt.show()
+    from Models.code_coreshellcyl import GpuCoreShellCylinder as gpumodel
+
+    plot(data, model, gpumodel, N, pars)
+
+def trellipse(N=1):
+
+    data, N = newlen(N)
+
+    pars = dict(scale=0.08, semi_axisA=15, semi_axisB=20, semi_axisC=500,
+                sldEll=7.105e-6, sldSolv=.291e-6,
+                background=5, axis_theta=0, axis_phi=0, axis_psi=0,
+                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
+
+    plot(data, model, gpumodel, N, pars)
+
+def lamellar(N=1):
+
+    data, N = newlen(N)
+
+    pars = dict(scale=0.08,
+                bi_thick=19.2946,
+                sld_bi=5.38e-6,sld_sol=7.105e-6,
+                background=0.003,
+                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
+
+    plot(data, model, gpumodel, N, pars)
+
+def cap(N=1):
+
+    data, N = newlen(N)
+
+    pars = dict(scale=.08, rad_cyl=20, rad_cap=40, len_cyl=400,
+                sld_capcyl=1e-6, sld_solv=6.3e-6,
+                background=0, 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
+
+    plot(data, model, gpumodel, N, pars)
+
 
 if __name__ == "__main__":
-    coreshell()
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+    cyl(2)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: fit.py
===================================================================
--- fit.py	(revision 8cdb9f16e9c711b409c9c462f8d7c27dddf0b001)
+++ fit.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -6,5 +6,6 @@
 from Models.code_capcyl import GpuCapCylinder
 from Models.code_coreshellcyl import GpuCoreShellCylinder
-from Models.code_cylinder import GpuCylinder, OneDGpuCylinder
+from Models.code_cylinder_f import GpuCylinder
+#from Models.code_cylinder import GpuCylinder, OneDGpuCylinder
 from Models.code_ellipse import GpuEllipse
 from Models.code_lamellar import GpuLamellar
@@ -13,9 +14,9 @@
 """ IMPORT THE DATA USED """
 #data = load_data('December/Tangential/Sector0/DEC07133.ABS')
-data = load_data('December/DEC07102.DAT')
+data = load_data('December/DEC07098.DAT')
 
 """ SET INNER BEAM STOP, OUTER RING, AND MASK HALF OF THE DATA """
-set_beam_stop(data, 0.00669)#, outer=0.025)
-set_top(data, -.018)
+set_beam_stop(data, 0.004)#, outer=0.025)
+#set_top(data, -.018)
 #set_half(data, 'right')
 
@@ -41,5 +42,5 @@
     radius_b_pd=.000128, radius_b_pd_n=1, radius_b_pd_nsigma=0,
     axis_phi_pd=2.63698e-05, axis_phi_pd_n=20, axis_phi_pd_nsigma=0,
-    dtype='float')
+    dtype='float32')
 
 
@@ -59,5 +60,5 @@
     background=0.003,
     bi_thick_pd= 0.37765, bi_thick_pd_n=10, bi_thick_pd_nsigma=3,
-    dtype='float')
+    dtype='float32')
 
     # SET THE FITTING PARAMETERS
@@ -70,19 +71,27 @@
 
 if 1:
-    model = SasModel(data, GpuCylinder,
-    scale=0.0104, radius=92.5, length=798.3,
-    sldCyl=.29e-6, sldSolv=7.105e-6, background=5,
-    cyl_theta=0, cyl_phi=0,
-    cyl_theta_pd=22.11, cyl_theta_pd_n=20, cyl_theta_pd_nsigma=3,
-    radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3,
-    length_pd=0.493, length_pd_n=10, length_pd_nsigma=3,
-    cyl_phi_pd=0, cyl_phi_pd_n=1, cyl_phi_pd_nsigma=3,
-    dtype='float')
+    """
+    pars = dict(scale=0.0023, radius=92.5, length=798.3,
+        sldCyl=.29e-6, sldSolv=7.105e-6, background=5,
+        cyl_theta=0, cyl_phi=0,
+        cyl_theta_pd=22.11, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3,
+        radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3,
+        length_pd=0.493, length_pd_n=10, length_pd_nsigma=3,
+        cyl_phi_pd=0, cyl_phi_pd_n=5, cyl_phi_pd_nsigma=3,)
+        """
+    pars = dict(
+        scale=.01, radius=64.1, length=66.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1,
+        cyl_theta=90, cyl_phi=0,
+        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
+        length_pd=0.1,length_pd_n=5, length_pd_nsigma=3,
+        cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3,
+        cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3)
+    model = SasModel(data, GpuCylinder, dtype="float32", **pars)
 
 
 
     # SET THE FITTING PARAMETERS
-    #model.radius.range(1, 500)
-    #model.length.range(1, 4000)
+    model.radius.range(1, 500)
+    model.length.range(1, 4000)
     #model.cyl_theta.range(-90,100)
     #model.cyl_theta_pd.range(0, 90)
@@ -106,5 +115,5 @@
     axis_theta_pd=1, axis_theta_pd_n=10, axis_theta_pd_nsigma=3,
     axis_phi_pd=0.1, axis_phi_pd_n=1, axis_phi_pd_nsigma=1,
-    dtype='float')
+    dtype='float32')
 
     # SET THE FITTING PARAMETERS
@@ -124,13 +133,13 @@
 if 0:
     model = SasModel(data, GpuCapCylinder,
-    scale=1, rad_cyl=20, rad_cap=40, length=400,
+    scale=.08, rad_cyl=20, rad_cap=40, len_cyl=400,
     sld_capcyl=1e-6, sld_solv=6.3e-6,
     background=0, theta=0, phi=0,
-    rad_cyl_pd=.1, rad_cyl_pd_n=1, rad_cyl_pd_nsigma=0,
-    rad_cap_pd=.1, rad_cap_pd_n=1, rad_cap_pd_nsigma=0,
-    length_pd=.1, length_pd_n=1, length_pd_nsigma=0,
+    rad_cyl_pd=.1, rad_cyl_pd_n=5, rad_cyl_pd_nsigma=3,
+    rad_cap_pd=.1, rad_cap_pd_n=5, rad_cap_pd_nsigma=3,
+    len_cyl_pd=.1, len_cyl_pd_n=1, len_cyl_pd_nsigma=0,
     theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0,
     phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0,
-    dtype='float')
+    dtype='float32')
 
     model.scale.range(0, 1)
@@ -147,5 +156,5 @@
     axisA_pd=.1, axisA_pd_n=1, axisA_pd_nsigma=0,
     axisB_pd=.1, axisB_pd_n=1, axisB_pd_nsigma=0,
-    axisC_pd=.1, axisC_pd_n=1, axisC_pd_nsigma=0, dtype='float')
+    axisC_pd=.1, axisC_pd_n=1, axisC_pd_nsigma=0, dtype='float32')
 
     # SET THE FITTING PARAMETERS
Index: sasmodel.py
===================================================================
--- sasmodel.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ sasmodel.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
@@ -2,4 +2,5 @@
 # -*- coding: utf-8 -*-
 
+import datetime
 import numpy as np
 import pyopencl as cl
@@ -8,4 +9,13 @@
 from sans.dataloader.manipulations import Ringcut, Boxcut
 
+
+TIC = None
+def tic():
+    global TIC
+    TIC = datetime.datetime.now()
+
+def toc():
+    now = datetime.datetime.now()
+    return (now-TIC).total_seconds()
 
 def load_data(filename):
@@ -21,4 +31,6 @@
     if np.dtype(dtype) == np.dtype('float32'):
         header = """\
+
+#define REAL(x) (x##f)
 #define real float
 """
@@ -26,4 +38,5 @@
         header = """\
 #pragma OPENCL EXTENSION cl_khr_fp64: enable
+#define REAL(x) (x)
 #define real double
 """
