Index: Kernel/Kernel-CoreShellCylinder_f.cpp
===================================================================
--- Kernel/Kernel-CoreShellCylinder_f.cpp	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
+++ Kernel/Kernel-CoreShellCylinder_f.cpp	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
@@ -0,0 +1,128 @@
+#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 subone, real subtwo, real radius, real length, real thickness, real weight, real axis_theta, real axis_phi);
+real f(real qx, real qy, real subone, real subtwo, real radius, real length, real thickness, real weight, real axis_theta, real axis_phi)
+{
+    const real q = sqrt(qx*qx+qy*qy);
+    const real theta = axis_theta*RADIANS;
+    real alpha = acos(cos(theta)*cos(axis_phi*RADIANS)*qx/q + sin(theta)*qy/q);
+
+    if (alpha == 0.0){
+    alpha = 1.0e-26;
+    }
+
+    real si1=0; real si2=0; real be1=0; real be2=0;
+
+    const real vol2 = M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness);
+
+    const real besarg1 = q*radius*sin(alpha);
+    const real besarg2 = q*(radius+thickness)*sin(alpha);
+    const real sinarg1 = q*length/2.0*cos(alpha);
+    const real sinarg2 = q*(length/2.0+thickness)*cos(alpha);
+
+    if (besarg1 == 0.0){be1 = 0.5;}
+    else{be1 = NR_BessJ1(besarg1)/besarg1;}
+
+    if (besarg2 == 0.0){be2 = 0.5;}
+    else{be2 = NR_BessJ1(besarg2)/besarg2;}
+
+    if (sinarg1 == 0.0){si1 = 1.0;}
+    else{si1 = sin(sinarg1)/sinarg1;}
+
+    if (sinarg2 == 0.0){si2 = 1.0;}
+    else{si2 = sin(sinarg2)/sinarg2;}
+
+    const real tt = 2.0*vol2*(subone)*si2*be2+2.0*(M_PI*radius*radius*(length))*(subtwo)*si1*be1;
+
+    real answer = (tt*tt)*sin(alpha)/fabs(sin(alpha));
+    //answer *= answer/(PI_E8*(radius+thickness)*(radius+thickness)*(length+2.0*thickness));
+
+    return weight*answer*answer/PI_E8;
+
+}
+
+kernel void CoreShellCylinderKernel(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 subone, const real subtwo,
+    const int Nradius, const int Nlength, const int Nthick, 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+Nthick+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];
+                for (int ti=0; ti < Nthick; ti++) {
+                    const real tv = loops[2*(ti+Nradius+Nlength)]; //////////////////////////
+                    const real tw = loops[2*(ti+Nradius+Nlength)+1];
+                    for (int thi=0; thi < Ntheta; thi++) {
+                        const real thv = loops[2*(thi+Nradius+Nlength+Nthick)];
+                        const real thw = loops[2*(thi+Nradius+Nlength+Nthick)+1];
+                        // #pragma unroll
+                        for (int phi=0; phi < Nphi; phi++) {
+                            const real phv = loops[2*(phi+Nradius+Nlength+Ntheta+Nthick)];
+                            const real phw = loops[2*(phi+Nradius+Nlength+Ntheta+Nthick)+1];
+
+                            const real weight = rw*lw*tw*thw*phw;
+                            //ret += qxi + qyi + sub + rv + lv + weight + thv + phv;
+                            if (weight > cutoff) {
+                                ret += f(qxi, qyi, subone, subtwo, rv, lv, tv, weight, thv, phv);
+                                norm += weight;
+                                vol += rw*lw*tw*(rv+tv)*(rv+tv)*(lv+2.0*tv);
+                                norm_vol += rw*lw*tw;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        //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: Kernel/Kernel-Cylinder_f.cpp
===================================================================
--- Kernel/Kernel-Cylinder_f.cpp	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ Kernel/Kernel-Cylinder_f.cpp	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
@@ -85,6 +85,4 @@
                 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)];
@@ -100,4 +98,6 @@
                             ret += f(qxi, qyi, sub, rv, lv, weight, thv, phv);
                             norm += weight;
+                            vol += rw*lw*rv*rv*lv;
+                            norm_vol += rw*lw;
                         }
                     }
@@ -105,5 +105,5 @@
             }
         }
-        if (Ntheta>1) norm = norm/(M_PI/2);
+        //if (Ntheta>1) norm = norm/(M_PI/2);
         if (vol != REAL(0.0) && norm_vol != REAL(0.0)) {
             ret *= norm_vol/vol;
Index: Kernel/Kernel-Ellipse_f.cpp
===================================================================
--- Kernel/Kernel-Ellipse_f.cpp	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
+++ Kernel/Kernel-Ellipse_f.cpp	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
@@ -73,6 +73,4 @@
                 const real rbv = loops[2*(bi+Nradius_a)];
                 const real rbw = loops[2*(bi+Nradius_a)+1];
-                vol += raw*rbw*rav*rbv*rbv*Ntheta*Nphi;
-                norm_vol += raw*rbw*Ntheta*Nphi;
                 for (int thi=0; thi < Ntheta; thi++) {
                     const real thv = loops[2*(thi+Nradius_a+Nradius_b)];
@@ -88,4 +86,6 @@
                             ret += f(qxi, qyi, sub, rav, rbv, weight, thv, phv);
                             norm += weight;
+                            vol += raw*rbw*rav*rbv*rbv;
+                            norm_vol += raw*rbw;
                         }
                     }
Index: Models/code_coreshellcyl_f.py
===================================================================
--- Models/code_coreshellcyl_f.py	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
+++ Models/code_coreshellcyl_f.py	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
@@ -0,0 +1,92 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+
+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 GpuCoreShellCylinder(object):
+    PARS = {'scale':1, 'radius':1, 'thickness':1, 'length':1, 'core_sld':1e-6, 'shell_sld':-1e-6, 'solvent_sld':0,
+            'background':0, 'axis_theta':0, 'axis_phi':0}
+    PD_PARS = ['radius', 'length', 'thickness', 'axis_phi', 'axis_theta']
+
+    def __init__(self, qx, qy, dtype='float32', cutoff=1e-5):
+        #create context, queue, and build program
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('Kernel/NR_BessJ1.cpp').read()+"\n"+open('Kernel/Kernel-CoreShellCylinder_f.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+        self.cutoff = cutoff
+
+        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, mf.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.CoreShellCylinderKernel(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['shell_sld']-pars['solvent_sld']),
+                        real(pars['core_sld']-pars['shell_sld']),
+                        *[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 CPUCoreShellCylinder(GpuCoreShellCylinder):
+    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/coreshellcylinder.so')
+        self.fn = self.dll['CoreShellCylinderKernel']
+        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_double,
+            c_int, 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['shell_sld']-pars['solvent_sld'], pars['core_sld']-pars['shell-sld'], *loop_lengths)
+
+        return self.res
+
+
+def make_loops(pars, dtype='double'):
+    # 0.2 ms on sparkle to form the final loops
+    radius, length, thickness, axis_theta, axis_phi = \
+       [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+         for base in GpuCoreShellCylinder.PD_PARS]
+    parts = [
+        radius.get_weights(pars['radius'], 0, 10000, True),
+        length.get_weights(pars['length'], 0, 10000, True),
+        thickness.get_weights(pars['thickness'], 0, 10000, True),
+        axis_theta.get_weights(pars['axis_theta'], -np.inf, np.inf, False),
+        axis_phi.get_weights(pars['axis_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]
Index: compare.py
===================================================================
--- compare.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ compare.py	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
@@ -117,5 +117,5 @@
 
     model = sasview_model('CoreShellCylinder', **pars)
-    from Models.code_coreshellcyl import GpuCoreShellCylinder as gpumodel
+    from Models.code_coreshellcyl_f import GpuCoreShellCylinder as gpumodel
 
     plot(data, model, gpumodel, N, pars)
@@ -176,32 +176,32 @@
 
 if __name__ == "__main__":
-    cyl(2)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+    coreshell(1)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: fit.py
===================================================================
--- fit.py	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
+++ fit.py	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
@@ -5,5 +5,5 @@
 from sasmodel import SasModel, load_data, set_beam_stop, set_half, set_top
 from Models.code_capcyl import GpuCapCylinder
-from Models.code_coreshellcyl import GpuCoreShellCylinder
+from Models.code_coreshellcyl_f import GpuCoreShellCylinder
 from Models.code_cylinder_f import GpuCylinder
 #from Models.code_cylinder import GpuCylinder, OneDGpuCylinder
@@ -18,5 +18,5 @@
 """ SET INNER BEAM STOP, OUTER RING, AND MASK HALF OF THE DATA """
 set_beam_stop(data, 0.004)#, outer=0.025)
-#set_top(data, -.018)
+set_top(data, -.018)
 #set_half(data, 'right')
 
@@ -31,5 +31,5 @@
 
 
-if 1:
+if 0:
     model = SasModel(data, GpuEllipse,
     scale=0.08,
@@ -104,9 +104,9 @@
 
 
-if 0:
+if 1:
     model = SasModel(data, GpuCoreShellCylinder,
-    scale= .00031, radius=19.5, thickness=30, length=22,
+    scale= .031, radius=19.5, thickness=30, length=22,
     core_sld=7.105e-6, shell_sld=.291e-6, solvent_sld=7.105e-6,
-    background=0.2, axis_theta=0, axis_phi=0,
+    background=0, axis_theta=0, axis_phi=0,
 
     radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3,
Index: sasmodel.py
===================================================================
--- sasmodel.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ sasmodel.py	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
@@ -101,4 +101,5 @@
     vmin = min(mdata.min(), mtheory.min())
     vmax = max(mdata.max(), mtheory.max())
+    print np.exp(np.mean(mtheory)), np.std(mtheory),np.max(mtheory),np.min(mtheory)
 
     plt.subplot(1, 3, 1)
