Index: Kernel/Kernel-Ellipse.cpp
===================================================================
--- Kernel/Kernel-Ellipse.cpp	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ Kernel/Kernel-Ellipse.cpp	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
@@ -1,4 +1,6 @@
-__kernel void EllipsoidKernel(const real radius_a_weight, const real radius_b_weight, const real axis_theta_weight,
-const real axis_phi_weight, const real scale, const real radius_a, const real radius_b, const real sub, const real axis_theta, const real axis_phi, __global const real *qx,
+__kernel void EllipsoidKernel(const real radius_a_weight, const real radius_b_weight,
+const real axis_theta_weight,
+const real axis_phi_weight, const real scale, const real radius_a, const real radius_b,
+ const real sub, const real axis_theta, const real axis_phi, __global const real *qx,
 __global const real *qy, __global real *_ptvalue, const int length, const int size)
 {
Index: Kernel/Kernel-Ellipse_f.cpp
===================================================================
--- Kernel/Kernel-Ellipse_f.cpp	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
+++ Kernel/Kernel-Ellipse_f.cpp	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
@@ -0,0 +1,101 @@
+#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"
+#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 radius_a, real radius_b, real scale, real axis_theta, real axis_phi);
+real f(real qx, real qy, real sub, real radius_a, real radius_b, real scale, real axis_theta, real axis_phi)
+{
+    real ret = 0;
+    const real q = sqrt(qx*qx + qy*qy);
+    const real theta = axis_theta*RADIANS;
+    const real cos_val = cos(theta)*cos(axis_phi*RADIANS)*(qx/q) + sin(theta)*(qy/q);
+
+    const real arg = q*radius_b*sqrt(1.0+(cos_val*cos_val*(((radius_a*radius_a/(radius_b*radius_b))-1.0))));
+    if(arg == 0.0){
+     ret = 1.0/3.0;
+    }
+    else{
+     ret = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
+    }
+    ret*=ret*9.0*sub*sub*4.0/3.0*acos(-1.0)*radius_b*radius_b*radius_a*(1.0e8);
+
+    return radius_a*scale*ret*pow(radius_b, 2);
+}
+
+kernel void EllipsoidKernel(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_a, const int Nradius_b, 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_a+Nradius_b+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 ai=0; ai < Nradius_a; ai++) {
+            const real rav = loops[2*ai];
+            const real raw = loops[2*ai+1];
+            for (int bi=0; bi < Nradius_b; bi++) {
+                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)];
+                    const real thw = loops[2*(thi+Nradius_a+Nradius_b)+1];
+                    // #pragma unroll
+                    for (int phi=0; phi < Nphi; phi++) {
+                        const real phv = loops[2*(phi+Nradius_a+Nradius_b+Ntheta)];
+                        const real phw = loops[2*(phi+Nradius_a+Nradius_b+Ntheta)+1];
+
+                        const real weight = raw*rbw*thw*phw;
+                        //ret += qxi + qyi + sub + rav + rbv + weight + thv + phv;
+                        if (weight > cutoff) {
+                            ret += f(qxi, qyi, sub, rav, rbv, 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_cylinder_f.py
===================================================================
--- Models/code_cylinder_f.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ Models/code_cylinder_f.py	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
@@ -2,5 +2,5 @@
 # -*- coding: utf-8 -*-
 
-import sys
+
 import ctypes
 from ctypes import c_int, c_double, c_void_p
@@ -45,4 +45,5 @@
         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,
@@ -51,4 +52,5 @@
             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]
@@ -85,5 +87,5 @@
         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),
+        phi.get_weights(pars['cyl_phi'], -np.inf, np.inf, False),
         ]
     # Make sure that weights are normalized to peaks at 1 so that
Index: Models/code_ellipse_f.py
===================================================================
--- Models/code_ellipse_f.py	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
+++ Models/code_ellipse_f.py	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
@@ -0,0 +1,89 @@
+#!/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 GpuEllipse(object):
+    PARS = {
+    'scale':1, 'radius_a':1, 'radius_b':1, 'sldEll':1e-6, 'sldSolv':0, 'background':0, 'axis_theta':0, 'axis_phi':0,
+    }
+    PD_PARS = ['radius_a', 'radius_b', 'axis_theta', 'axis_phi']
+
+    def __init__(self, qx, qy, dtype='float32', cutoff=1e-5):
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('Kernel/Kernel-Ellipse_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.EllipsoidKernel(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['sldEll']-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(GpuEllipse):
+    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/ellipsoid.so')
+        self.fn = self.dll['EllipsoidKernel']
+        self.fn.argtypes = [
+            c_void_p, c_void_p, c_void_p, c_int,
+            c_void_p, 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['sldEll']-pars['sldSolv'], *loop_lengths)
+
+        return self.res
+
+def make_loops(pars, dtype='double'):
+    # 0.2 ms on sparkle to form the final loops
+    radius_a, radius_b, theta, phi = \
+       [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+         for base in GpuEllipse.PD_PARS]
+    parts = [
+        radius_a.get_weights(pars['radius_a'], 0, 10000, True),
+        radius_b.get_weights(pars['radius_b'], 0, 10000, True),
+        theta.get_weights(pars['axis_theta'], -np.inf, np.inf, False),
+        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: fit.py
===================================================================
--- fit.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ fit.py	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
@@ -8,5 +8,5 @@
 from Models.code_cylinder_f import GpuCylinder
 #from Models.code_cylinder import GpuCylinder, OneDGpuCylinder
-from Models.code_ellipse import GpuEllipse
+from Models.code_ellipse_f import GpuEllipse
 from Models.code_lamellar import GpuLamellar
 from Models.code_triaxialellipse import GpuTriEllipse
@@ -31,5 +31,5 @@
 
 
-if 0:
+if 1:
     model = SasModel(data, GpuEllipse,
     scale=0.08,
@@ -41,14 +41,14 @@
     radius_a_pd=0.222296, radius_a_pd_n=1, radius_a_pd_nsigma=0,
     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,
+    axis_phi_pd=2.63698e-05, axis_phi_pd_n=20, axis_phi_pd_nsigma=3,
     dtype='float32')
 
 
     # SET THE FITTING PARAMETERS
-    model.radius_a.range(15, 1000)
-    model.radius_b.range(15, 1000)
+    #model.radius_a.range(15, 1000)
+    #model.radius_b.range(15, 1000)
     #model.axis_theta_pd.range(0, 360)
     #model.background.range(0,1000)
-    #model.scale.range(0, 1)
+    model.scale.range(0, 1)
 
 
@@ -70,5 +70,5 @@
 
 
-if 1:
+if 0:
     """
     pars = dict(scale=0.0023, radius=92.5, length=798.3,
