Index: Capcyl_Gauss.py
===================================================================
--- Capcyl_Gauss.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
+++ Capcyl_Gauss.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -0,0 +1,98 @@
+Gauss76Z = {.999505948362153*(-1.0),.997397786355355*(-1.0),.993608772723527*(-1.0),.988144453359837*(-1.0),
+               .981013938975656*(-1.0),.972229228520377*(-1.0),.961805126758768*(-1.0),.949759207710896*(-1.0),
+               .936111781934811*(-1.0),.92088586125215*(-1.0),.904107119545567*(-1.0),.885803849292083*(-1.0),
+               .866006913771982*(-1.0),.844749694983342*(-1.0),.822068037328975*(-1.0),.7980001871612*(-1.0),
+               .77258672828181*(-1.0),.74587051350361*(-1.0),.717896592387704*(-1.0),.688712135277641*(-1.0),
+               .658366353758143*(-1.0),.626910417672267*(-1.0),.594397368836793*(-1.0),.560882031601237*(-1.0),
+               .526420920401243*(-1.0),.491072144462194*(-1.0),.454895309813726*(-1.0),.417951418780327*(-1.0),
+               .380302767117504*(-1.0),.342012838966962*(-1.0),.303146199807908*(-1.0),.263768387584994*(-1.0),
+               .223945802196474*(-1.0),.183745593528914*(-1.0),.143235548227268*(-1.0),.102483975391227*(-1.0),
+               .0615595913906112*(-1.0),.0205314039939986*(-1.0),-.0205314039939986*(-1.0),-.0615595913906112*(-1.0),
+               -.102483975391227*(-1.0),-.143235548227268*(-1.0),-.183745593528914*(-1.0),-.223945802196474*(-1.0),
+               -.263768387584994*(-1.0),-.303146199807908*(-1.0),-.342012838966962*(-1.0),-.380302767117504*(-1.0),
+               -.417951418780327*(-1.0),-.454895309813726*(-1.0),-.491072144462194*(-1.0),-.526420920401243*(-1.0),
+               -.560882031601237*(-1.0),-.594397368836793*(-1.0),-.626910417672267*(-1.0),-.658366353758143*(-1.0),
+               -.688712135277641*(-1.0),-.717896592387704*(-1.0),-.74587051350361*(-1.0),-.77258672828181*(-1.0),
+               -.7980001871612*(-1.0),	-.822068037328975*(-1.0),-.844749694983342*(-1.0),-.866006913771982*(-1.0),
+               -.885803849292083*(-1.0),-.904107119545567*(-1.0),-.92088586125215*(-1.0),-.936111781934811*(-1.0),
+               -.949759207710896*(-1.0),-.961805126758768*(-1.0),-.972229228520377*(-1.0),-.981013938975656*(-1.0),
+               -.988144453359837*(-1.0),-.993608772723527*(-1.0),-.997397786355355*(-1.0),-.999505948362153*(-1.0) }
+
+Gauss76Wt= {
+.00126779163408536,
+.00294910295364247,
+.00462793522803742,
+.00629918049732845,
+.00795984747723973,
+.00960710541471375,
+.0112381685696677,
+.0128502838475101,
+.0144407317482767,
+.0160068299122486,
+.0175459372914742,
+.0190554584671906,
+.020532847967908,
+.0219756145344162,
+.0233813253070112,
+.0247476099206597,
+.026072164497986,
+.0273527555318275,
+.028587223650054,
+.029773487255905,
+.0309095460374916,
+.0319934843404216,
+.0330234743977917,
+.0339977794120564,
+.0349147564835508,
+.0357728593807139,
+.0365706411473296,
+.0373067565423816,
+.0379799643084053,
+.0385891292645067,
+.0391332242205184,
+.0396113317090621,
+.0400226455325968,
+.040366472122844,
+.0406422317102947,
+.0408494593018285,
+.040987805464794,
+.0410570369162294,
+.0410570369162294,
+.040987805464794,
+.0408494593018285,
+.0406422317102947,
+.040366472122844,
+.0400226455325968,
+.0396113317090621,
+.0391332242205184,
+.0385891292645067,
+.0379799643084053,
+.0373067565423816,
+.0365706411473296,
+.0357728593807139,
+.0349147564835508,
+.0339977794120564,
+.0330234743977917,
+.0319934843404216,
+.0309095460374916,
+.029773487255905,
+.028587223650054,
+.0273527555318275,
+.026072164497986,
+.0247476099206597,
+.0233813253070112,
+.0219756145344162,
+.020532847967908,
+.0190554584671906,
+.0175459372914742,
+.0160068299122486,
+.0144407317482767,
+.0128502838475101,
+.0112381685696677,
+.00960710541471375,
+.00795984747723973,
+.00629918049732845,
+.00462793522803742,
+.00294910295364247,
+.00126779163408536
+}
Index: Capcyl_Kfun.cpp
===================================================================
--- Capcyl_Kfun.cpp	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
+++ Capcyl_Kfun.cpp	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -14,5 +14,6 @@
 	return(val);
 }
-const real Gauss76Z[76]={
+/*
+real Gauss76Z[76]={
 	 .999505948362153*(-1.0),		//0
 	 .997397786355355*(-1.0),
@@ -92,5 +93,5 @@
 	-.999505948362153*(-1.0)		//75
 };
-const real Gauss76Wt[76]={
+real Gauss76Wt[76]={
 	.00126779163408536,		//0
 	.00294910295364247,
@@ -170,2 +171,3 @@
 	.00126779163408536		//75 (indexed from 0)
 };
+*/
Index: Kernel-CapCyl.cpp
===================================================================
--- Kernel-CapCyl.cpp	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
+++ Kernel-CapCyl.cpp	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -1,6 +1,6 @@
-__kernel void CapCylinderKernel(__global const float *qx, __global const float *qy, __global float *__ptvalue, __global float *vol_i
+__kernel void CapCylinderKernel(__global const float *qx, __global const float *qy, __global float *_ptvalue, __global float *vol_i,
 const float rad_cyl, const float rad_cap, const float length, const float thet, const float ph, const float sub,
 const float scale, const float phi_weight, const float theta_float, const float rad_cap_weight, const float rad_cyl_weight,
-const float length_weight, const int total, const int size)
+const float length_weight, const int total, const int size, __const float Gauss76Wt, __const float Gauss76Z)
 //ph is phi, sub is sldc-slds, thet is theta
 {
@@ -8,5 +8,5 @@
     if(i < total)
     {
-        float q = sqrt(qx[i]*qx[i] + qy[i]*qy[i])
+        float q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
         float pi = 4.0*atan(1.0);
         float theta = thet*pi/180.0;
@@ -16,11 +16,11 @@
         float cos_val = cyl_x*qx[i]/q + cyl_y*qy[i]/q;
         float alpha = acos(cos_val);
-        float yyy=0; float ans1=0; float ans2=0; float y=0; float xx=0; float ans=0; float zij=0; float be=0
+        float yyy=0; float ans1=0; float ans2=0; float y=0; float xx=0; float ans=0; float zij=0; float be=0; float summj=0;
 
         float hDist = -1.0*sqrt(fabs(rad_cap*rad_cap-rad_cyl*rad_cyl));
-        vol_i[i] = pi*rad_cyl*rad_cyl*len_cyl+2.0*pi/3.0*((rad_cap-hDist)*(rad_cap-hDist)*(2*rad_cap+hDist));
+        vol_i[i] = pi*rad_cyl*rad_cyl*length+2.0*pi/3.0*((rad_cap-hDist)*(rad_cap-hDist)*(2*rad_cap+hDist));
         float vaj = -1.0*hDist/rad_cap;
 
-        for(j=0;j<76;j++) //the 76 corresponds to the Gauss constants
+        for(int j=0;j<76;j++) //the 76 corresponds to the Gauss constants
         {
             zij = (Gauss76Z[j]*(1.0-vaj)+vaj+1.0)/2.0;
@@ -45,5 +45,5 @@
         }
         float answer=yyy*yyy*1.0e8*sub*sub*scale/pi*rad_cyl*rad_cyl*length+2.0*pi*(2.0*rad_cap*rad_cap*rad_cap/3.0+rad_cap*rad_cap*hDist-hDist*hDist*hDist/3.0);
-        answer/=sin(alpha)
+        answer/=sin(alpha);
 
         _ptvalue[i] = rad_cyl_weight*length_weight*rad_cap_weight*theta_weight*phi_weight*vol_i[i]*answer
Index: Kernel-CoreShellCylinder.cpp
===================================================================
--- Kernel-CoreShellCylinder.cpp	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
+++ Kernel-CoreShellCylinder.cpp	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -11,9 +11,7 @@
         float pi = 4.0*atan(1.0);
         float theta = axis_theta*pi/180.0;
-        float phi = axis_phi*pi/180.0;
-        float cyl_x = cos(theta)*cos(phi);
+        float cyl_x = cos(theta)*cos(axis_phi*pi/180.0);
         float cyl_y = sin(theta);
-        float cos_val = cyl_x*qx[i]/q + cyl_y*qx[i]/q;
-        float alpha = acos(cos_val);
+        float alpha = acos(cyl_x*qx[i]/q + cyl_y*qx[i]/q);
 
         if (alpha == 0.0){
@@ -23,7 +21,4 @@
         float si1=0; float si2=0; float be1=0; float be2=0;
 
-        float dr1 = core_sld-shell_sld;
-        float dr2 = shell_sld-solvent_sld;
-        float vol1 = pi*radius*radius*(length);
         float vol2 = pi*(radius+thickness)*(radius+thickness)*(length+2.0*thickness);
 
@@ -54,9 +49,7 @@
             si2 = sin(sinarg2)/sinarg2;
         }
+        float tt = 2.0*vol2*(shell_sld-solvent_sld)*si2*be2 + 2.0*(pi*radius*radius*(length))*(core_sld-shell_sld)*si1*be1;
 
-        float t1 = 2.0*vol1*dr1*si1*be1;
-        float t2 = 2.0*vol2*dr2*si2*be2;
-
-        float answer = ((t1+t2)*(t1+t2))*sin(alpha)/fabs(sin(alpha));
+        float answer = (tt*tt)*sin(alpha)/fabs(sin(alpha));
         float vol=pi*(radius+thickness)*(radius+thickness)*(length+2.0*thickness);
         answer = answer/vol*1.0e8*scale;
Index: Kernel-Lamellar.cpp
===================================================================
--- Kernel-Lamellar.cpp	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ Kernel-Lamellar.cpp	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -1,8 +1,5 @@
-#ifndef real
-# define real float
-#endif
 
 __kernel void LamellarKernel(__global const real *qx, global const real *qy, __global real *ret, const real bi_thick,
- const real scale, const real sub, const real background, const int length)
+ const real scale, const real sub, const int length)
 {
     int i = get_global_id(0);
@@ -12,6 +9,5 @@
         real pi = 4.0*atan(1.0);
         real Pq = 2.0*sub*(sub/q)/q*(1.0-cos(q*bi_thick));
-        ret[i] = 2.0*pi*scale*Pq/(q*q)/bi_thick*1.0e8;
-        ret[i] += background;
+        ret[i] = 2.0*pi*scale*Pq/(q*q)/bi_thick*1.0e8 ;
     }
 }
Index: Kernel-TriaxialEllipse.cpp
===================================================================
--- Kernel-TriaxialEllipse.cpp	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
+++ Kernel-TriaxialEllipse.cpp	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -0,0 +1,45 @@
+__kernel void TriaxialEllipseKernel(__global const real *qx, __global const real *qy, __global real *_ptvalue,
+const real sub, const real scale, const real axisA, const real axisB, const real axisC, const real axis_phi,
+const real axis_theta, const real axis_psi, const real axisA_weight, const real axisB_weight, const real axisC_weight,
+const real psi_weight, const real phi_weight, const real theta_weight, const int total, const int size)
+{
+    int i = get_global_id(0);
+    if(i < total)
+    {
+        real q = sqrt(qx[i]*qx[i]+qy[i]*qy[i]);
+        real q_x = qx[i]/q;
+        real q_y = qy[i]/q;
+        real pi = 4.0*atan(1.0);
+        
+        //convert angle degree to radian
+        real theta = axis_theta*pi/180.0;
+        real phi = axis_phi*pi/180.0;
+        real psi = axis_psi*pi/180.0;
+        real cos_val = cos(theta)*cos(phi)*q_x + sin(theta)*q_y;
+        
+        // calculate the axis of the ellipse wrt q-coord.
+        real cos_nu = (-cos(phi)*sin(psi)*sin(theta)+sin(phi)*cos(psi))*q_x + sin(psi)*cos(theta)*q_y;
+        real cos_mu = (-sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi))*q_x + cos(theta)*cos(psi)*q_y;
+
+        real answer=0;
+        real t = q*sqrt(axisA*axisA*cos_nu*cos_nu+axisB*axisB*cos_mu*cos_mu+axisC*axisC*cos_val*cos_val);
+
+        if (t==0.0){
+            answer = 1.0;
+        }
+        else{
+            answer  = 3.0*(sin(t)-t*cos(t))/(t*t*t);
+        }
+        answer*=answer*sub*sub*4.0*pi/3.0*axisA*axisB*axisC*1.0e8*scale;
+
+        _ptvalue[i] = axisA_weight*axisB_weight*axisC_weight*theta_weight*phi_weight*psi_weight*answer*axisA*axisB*axisC;
+        if (size>1)
+        {
+            _ptvalue[i] *= fabs(cos(theta*pi/180.0));
+        }
+    }
+}
+
+
+
+
Index: pcylcope.py
===================================================================
--- capcylcope.py	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
+++ 	(revision )
@@ -1,131 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-from math import asin
-import pyopencl as cl
-from weights import GaussianDispersion
-from sasmodel import card
-
-def set_precision(src, qx, qy, dtype):
-    qx = np.ascontiguousarray(qx, dtype=dtype)
-    qy = np.ascontiguousarray(qy, dtype=dtype)
-    if np.dtype(dtype) == np.dtype('float32'):
-        header = """\
-#define real float
-"""
-    else:
-        header = """\
-#pragma OPENCL EXTENSION cl_khr_fp64: enable
-#define real double
-"""
-    return header+src, qx, qy
-
-class GpuCapCylinder(object):
-    PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'length':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0,
-             'theta':0, 'phi':0}
-
-    PD_PARS = ['rad_cyl', 'length', 'rad_cap', 'theta', 'phi']
-
-    def __init__(self, qx, qy, dtype='float32'):
-
-        #create context, queue, and build program
-        ctx,_queue = card()
-        trala = open('NR_BessJ1.cpp').read()+"\n"+open('Capcyl_Kfun.cpp').read()+"\n"+open('Kernel-Cylinder.cpp').read()
-        src, qx, qy = set_precision(trala, qx, qy, dtype=dtype)
-
-        self.prg = cl.Program(ctx, open('Kernel-CapCyl.cpp').read()).build()
-
-        #buffers
-        mf = cl.mem_flags
-        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.WRITE_ONLY, qx.nbytes)
-        self.res = np.empty_like(self.qx)
-        self.vol_i = float(0.0)
-        self.vol_b = cl.Buffer(ctx, mf.WRITE_ONLY, self.vol_i.nbytes)
-
-    def eval(self, pars):
-
-        _ctx,queue = card()
-        rad_cyl,length,rad_cap,theta,phi = \
-            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
-             for base in GpuCapCylinder.PD_PARS]
-
-        rad_cyl.value, rad_cyl.weight = rad_cyl.get_weights(pars['rad_cyl'], 0, 1000, True)
-        rad_cap.value, rad_cap.weight = rad_cap.get_weights(pars['rad_cap'], 0, 1000, True)
-        length.value, length.weight = length.get_weights(pars['length'], 0, 1000, 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)
-
-        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
-        size = len(theta.weight)
-        sub = pars['sld_capcyl']-np.float32(['sld_solv'])
-
-        for i in xrange(len(rad_cyl.weight)):
-            for m in xrange(len(rad_cap.weight)):
-                for j in xrange(len(length.weight)):
-                    for k in xrange(len(theta.weight)):
-                        for l in xrange(len(phi.weight)):
-
-                            self.prg.CapCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                        self.vol_b, np.float32(rad_cyl.value[i]), np.float32(rad_cap.value[m]), np.float32(length.value[j]),
-                                        np.float32(theta.value[k]), np.float32(phi.value[l]), np.float32(sub), np.float32(pars['scale']),
-                                        np.float32(phi.weight[l]), np.float32(theta.weight[k]), np.float32(rad_cap.weight[m]),
-                                        np.float32(rad_cyl.weight[i]), np.float32(length.weight[j]), np.uint32(self.qx.size), np.uint32(size))
-
-                            cl.enqueue_copy(queue, self.res, self.res_b)
-                            cl.enqueue_copy(queue, self.vol_i, self.vol_b)
-
-                            sum += self.res
-                            vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*self.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]
-
-        if size > 1:
-            norm /= asin(1.0)
-
-        if vol != 0.0 and norm_vol != 0.0:
-            sum *= norm_vol/vol
-
-        return sum/norm + pars['background']
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: code_capcyl.py
===================================================================
--- code_capcyl.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
+++ code_capcyl.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -0,0 +1,137 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+from math import asin
+import pyopencl as cl
+from weights import GaussianDispersion
+from sasmodel import card
+from Capcyl_Gauss import Gauss76Wt, Gauss76Z
+
+def set_precision(src, qx, qy, dtype):
+    qx = np.ascontiguousarray(qx, dtype=dtype)
+    qy = np.ascontiguousarray(qy, dtype=dtype)
+    if np.dtype(dtype) == np.dtype('float32'):
+        header = """\
+#define real float
+"""
+    else:
+        header = """\
+#pragma OPENCL EXTENSION cl_khr_fp64: enable
+#define real double
+"""
+    return header+src, qx, qy
+
+class GpuCapCylinder(object):
+    PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'length':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0,
+             'theta':0, 'phi':0}
+
+    PD_PARS = ['rad_cyl', 'length', 'rad_cap', 'theta', 'phi']
+
+    def __init__(self, qx, qy, dtype='float32'):
+
+        #create context, queue, and build program
+        ctx,_queue = card()
+        trala = open('NR_BessJ1.cpp').read()+"\n"+open('Capcyl_Kfun.cpp').read()+"\n"+open('Kernel-Cylinder.cpp').read()
+        src, qx, qy = set_precision(trala, qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+
+        #buffers
+        mf = cl.mem_flags
+        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)
+        G, Z = Gauss76Wt, Gauss76Z
+        self.Gauss76W_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=G)
+        self.Gauss76Z_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=Z)
+        self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, qx.nbytes)
+        self.res = np.empty_like(self.qx)
+        self.vol_i = float(0.0)
+        self.vol_b = cl.Buffer(ctx, mf.WRITE_ONLY, self.vol_i.nbytes)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        rad_cyl,length,rad_cap,theta,phi = \
+            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+             for base in GpuCapCylinder.PD_PARS]
+
+        rad_cyl.value, rad_cyl.weight = rad_cyl.get_weights(pars['rad_cyl'], 0, 1000, True)
+        rad_cap.value, rad_cap.weight = rad_cap.get_weights(pars['rad_cap'], 0, 1000, True)
+        length.value, length.weight = length.get_weights(pars['length'], 0, 1000, 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)
+
+        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
+        size = len(theta.weight)
+        sub = pars['sld_capcyl']-np.float32(['sld_solv'])
+
+        for i in xrange(len(rad_cyl.weight)):
+            for m in xrange(len(rad_cap.weight)):
+                for j in xrange(len(length.weight)):
+                    for k in xrange(len(theta.weight)):
+                        for l in xrange(len(phi.weight)):
+
+                            self.prg.CapCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
+                                        self.vol_b, np.float32(rad_cyl.value[i]), np.float32(rad_cap.value[m]), np.float32(length.value[j]),
+                                        np.float32(theta.value[k]), np.float32(phi.value[l]), np.float32(sub), np.float32(pars['scale']),
+                                        np.float32(phi.weight[l]), np.float32(theta.weight[k]), np.float32(rad_cap.weight[m]),
+                                        np.float32(rad_cyl.weight[i]), np.float32(length.weight[j]), np.uint32(self.qx.size), np.uint32(size),
+                                        self.Gauss76W_b, self.Gauss76Z_b)
+
+                            cl.enqueue_copy(queue, self.res, self.res_b)
+                            cl.enqueue_copy(queue, self.vol_i, self.vol_b)
+
+                            sum += self.res
+                            vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*self.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]
+
+        if size > 1:
+            norm /= asin(1.0)
+
+        if vol != 0.0 and norm_vol != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm + pars['background']
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: code_coreshellcyl.py
===================================================================
--- code_coreshellcyl.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
+++ code_coreshellcyl.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -0,0 +1,86 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import math
+import pyopencl as cl
+from weights import GaussianDispersion
+from sasmodel import card
+
+def set_precision(src, qx, qy, dtype):
+    qx = np.ascontiguousarray(qx, dtype=dtype)
+    qy = np.ascontiguousarray(qy, dtype=dtype)
+    if np.dtype(dtype) == np.dtype('float32'):
+        header = """\
+#define real float
+"""
+    else:
+        header = """\
+#pragma OPENCL EXTENSION cl_khr_fp64: enable
+#define real double
+"""
+    return header+src, qx, qy
+
+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'):
+        #create context, queue, and build program
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('NR_BessJ1.cpp').read()+"\n"+open('Kernel-CoreShellCylinder.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+
+        #buffers
+        mf = cl.mem_flags
+        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.WRITE_ONLY, qx.nbytes)
+        self.res = np.empty_like(qx)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        radius, length, thickness, axis_phi, axis_theta = [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+                                     for base in GpuCoreShellCylinder.PD_PARS]
+
+        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 1000, True)
+        length.value, length.weight = length.get_weights(pars['length'], 0, 1000, True)
+        thickness.value, thickness.weight = thickness.get_weights(pars['thickness'], 0, 1000, True)
+        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False)
+        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False)
+
+        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
+        size = len(axis_theta.weight)
+
+        for i in xrange(len(radius.weight)):
+            for j in xrange(len(length.weight)):
+                for k in xrange(len(axis_theta.weight)):
+                    for l in xrange(len(axis_phi.weight)):
+                        for f in xrange(len(thickness.weight)):
+
+                            self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
+                                    np.float32(axis_theta.value[k]), np.float32(axis_phi.value[l]), np.float32(thickness.value[f]),
+                                    np.float32(length.value[j]), np.float32(radius.value[i]), np.float32(pars['scale']),
+                                    np.float32(radius.weight[i]), np.float32(length.weight[j]), np.float32(thickness.weight[f]),
+                                    np.float32(axis_theta.weight[k]), np.float32(axis_phi.weight[l]), np.float32(pars['core_sld']),
+                                    np.float32(pars['shell_sld']), np.float32(pars['solvent_sld']),np.uint32(size),
+                                    np.uint32(self.qx.size))
+                            cl.enqueue_copy(queue, self.res, self.res_b)
+
+                            sum += self.res
+                            vol += radius.weight[i]*length.weight[j]*thickness.weight[f]*pow(radius.value[i]+thickness.value[f],2)\
+                                   *(length.value[j]+2.0*thickness.value[f])
+                            norm_vol += radius.weight[i]*length.weight[j]*thickness.weight[k]
+                            norm += radius.weight[i]*length.weight[j]*thickness.weight[f]*axis_theta.weight[k]\
+                                    *axis_phi.weight[l]
+
+        if size>1:
+            norm /= math.asin(1.0)
+        if vol != 0.0 and norm_vol != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm + pars['background']
Index: code_cylinder.py
===================================================================
--- code_cylinder.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
+++ code_cylinder.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -0,0 +1,148 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import math
+import pyopencl as cl
+from weights import GaussianDispersion
+from sasmodel import card
+
+def set_precision(src, qx, qy, dtype):
+    qx = np.ascontiguousarray(qx, dtype=dtype)
+    qy = np.ascontiguousarray(qy, dtype=dtype)
+    if np.dtype(dtype) == np.dtype('float32'):
+        header = """\
+#define real float
+"""
+    else:
+        header = """\
+#pragma OPENCL EXTENSION cl_khr_fp64: enable
+#define real double
+"""
+    return header+src, qx, qy
+
+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'):
+
+        #create context, queue, and build program
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('NR_BessJ1.cpp').read()+"\n"+open('Kernel-Cylinder.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+        #buffers
+        mf = cl.mem_flags
+        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.WRITE_ONLY, qx.nbytes)
+        self.res = np.empty_like(self.qx)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        radius, length, cyl_theta, cyl_phi = \
+            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+             for base in GpuCylinder.PD_PARS]
+
+        #Get the weights for each
+        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 1000, True)
+        length.value, length.weight = length.get_weights(pars['length'], 0, 1000, True)
+        cyl_theta.value, cyl_theta.weight = cyl_theta.get_weights(pars['cyl_theta'], -90, 180, False)
+        cyl_phi.value, cyl_phi.weight = cyl_phi.get_weights(pars['cyl_phi'], -90, 180, False)
+
+        #Perform the computation, with all weight points
+        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
+        size = len(cyl_theta.weight)
+        sub = pars['sldCyl'] - pars['sldSolv']
+
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
+        #Loop over radius, length, theta, phi weight points
+        for i in xrange(len(radius.weight)):
+            for j in xrange(len(length.weight)):
+                for k in xrange(len(cyl_theta.weight)):
+                    for l in xrange(len(cyl_phi.weight)):
+
+
+                        self.prg.CylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(sub),
+                                           real(radius.value[i]), real(length.value[j]), real(pars['scale']),
+                                           real(radius.weight[i]), real(length.weight[j]), real(cyl_theta.weight[k]),
+                                           real(cyl_phi.weight[l]), real(cyl_theta.value[k]), real(cyl_phi.value[l]),
+                                           np.uint32(self.qx.size), np.uint32(size))
+                        cl.enqueue_copy(queue, self.res, self.res_b)
+                        sum += self.res
+                        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:
+            norm /= math.asin(1.0)
+        if vol != 0.0 and norm_vol != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm+pars['background']
+
+def demo():
+    from time import time
+    import matplotlib.pyplot as plt
+
+    #create qx and qy evenly spaces
+    qx = np.linspace(-.02, .02, 128)
+    qy = np.linspace(-.02, .02, 128)
+    qx, qy = np.meshgrid(qx, qy)
+
+    #saved shape of qx
+    r_shape = qx.shape
+
+    #reshape for calculation; resize as float32
+    qx = qx.flatten()
+    qy = qy.flatten()
+
+    pars = CylinderParameters(scale=1, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=0,
+                              cyl_theta=0, cyl_phi=0)
+    t = time()
+    result = GpuCylinder(qx, qy)
+    result.x = result.cylinder_fit(pars, r_n=10, t_n=10, l_n=10, p_n=10, r_w=.1, t_w=.1, l_w=.1, p_w=.1, sigma=3)
+    result.x = np.reshape(result.x, r_shape)
+    tt = time()
+
+    print("Time taken: %f" % (tt - t))
+
+    plt.pcolormesh(result.x)
+    plt.show()
+
+
+if __name__=="__main__":
+    demo()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: code_ellipse.py
===================================================================
--- code_ellipse.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
+++ code_ellipse.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -0,0 +1,113 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import math
+import pyopencl as cl
+from weights import GaussianDispersion
+
+
+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):
+
+        self.qx = np.asarray(qx, np.float32)
+        self.qy = np.asarray(qy, np.float32)
+        #create context, queue, and build program
+        self.ctx = cl.create_some_context()
+        self.queue = cl.CommandQueue(self.ctx)
+        self.prg = cl.Program(self.ctx, open('Kernel-Ellipse.cpp').read()).build()
+
+        #buffers
+        mf = cl.mem_flags
+        self.qx_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
+        self.qy_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
+        self.res_b = cl.Buffer(self.ctx, mf.WRITE_ONLY, qx.nbytes)
+        self.res = np.empty_like(self.qx)
+
+    def eval(self, pars):
+    #b_n = radius_b # want, a_n = radius_a # want, etc
+        radius_a, radius_b, axis_theta, axis_phi = \
+            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+             for base in GpuEllipse.PD_PARS]
+
+        radius_a.value, radius_a.weight = radius_a.get_weights(pars['radius_a'], 0, 1000, True)
+        radius_b.value, radius_b.weight = radius_b.get_weights(pars['radius_b'], 0, 1000, 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)
+
+        #Perform the computation, with all weight points
+        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
+        size = len(axis_theta.weight)
+        sub =  pars['sldEll'] - pars['sldSolv']
+
+        #Loop over radius weight points
+        for i in xrange(len(radius_a.weight)):
+            #Loop over length weight points
+            for j in xrange(len(radius_b.weight)):
+                #Average over theta distribution
+                for k in xrange(len(axis_theta.weight)):
+                    #Average over phi distribution
+                    for l in xrange(len(axis_phi.weight)):
+                        #call the kernel
+                        self.prg.EllipsoidKernel(self.queue, self.qx.shape, None, np.float32(radius_a.weight[i]),
+                                        np.float32(radius_b.weight[j]), np.float32(axis_theta.weight[k]),
+                                        np.float32(axis_phi.weight[l]), np.float32(pars['scale']), np.float32(radius_a.value[i]),
+                                        np.float32(radius_b.value[j]), np.float32(sub),np.float32(pars['background']),
+                                        np.float32(axis_theta.value[k]), np.float32(axis_phi.value[l]), self.qx_b, self.qy_b,
+                                        self.res_b, np.uint32(self.qx.size), np.uint32(len(axis_theta.weight)))
+                        #copy result back from buffer
+                        cl.enqueue_copy(self.queue, self.res, self.res_b)
+                        sum += self.res
+                        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
+        # integration (see documentation).
+        if size > 1:
+            norm /= math.asin(1.0)
+        if vol.any() != 0.0 and norm_vol.any() != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm+pars['background']
+
+
+def demo():
+    from time import time
+    import matplotlib.pyplot as plt
+
+    #create qx and qy evenly spaces
+    qx = np.linspace(-.02, .02, 128)
+    qy = np.linspace(-.02, .02, 128)
+    qx, qy = np.meshgrid(qx, qy)
+
+    #saved shape of qx
+    r_shape = qx.shape
+    #reshape for calculation; resize as float32
+    qx = qx.flatten()
+    qy = qy.flatten()
+
+    #int main
+    pars = EllipsoidParameters(.027, 60, 180, .297e-6, 5.773e-06, 4.9, 0, 90)
+
+    t = time()
+    result = GpuEllipse(qx, qy)
+    result.x = result.ellipsoid_fit(qx, qy, pars, b_n=35, t_n=35, a_n=1, p_n=1, sigma=3, b_w=.1, t_w=.1, a_w=.1, p_w=.1)
+    result.x = np.reshape(result.x, r_shape)
+    tt = time()
+    print("Time taken: %f" % (tt - t))
+
+    plt.pcolormesh(result.x)
+    plt.show()
+
+
+if __name__ == "__main__":
+    demo()
+
+
+
+
Index: code_lamellar.py
===================================================================
--- code_lamellar.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
+++ code_lamellar.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -0,0 +1,131 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import pyopencl as cl
+from weights import GaussianDispersion
+
+def set_precision(src, qx, qy, dtype):
+    qx = np.ascontiguousarray(qx, dtype=dtype)
+    qy = np.ascontiguousarray(qy, dtype=dtype)
+    if dtype == 'double':
+        header = """\
+#define real float
+"""
+    else:
+        header = """\
+#pragma OPENCL EXTENSION cl_khr_fp64: enable
+#define real double
+"""
+    return header+src, qx, qy
+
+
+class GpuLamellar(object):
+    PARS = {
+        'scale':1, 'bi_thick':1, 'sld_bi':1e-6, 'sld_sol':0, 'background':0,
+    }
+    PD_PARS = {'bi_thick'}
+    def __init__(self, qx, qy, dtype='float32'):
+
+        #create context, queue, and build program
+        self.ctx = cl.create_some_context()
+        self.queue = cl.CommandQueue(self.ctx)
+        src,qx,qy = set_precision(open('Kernel-Lamellar.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(self.ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+        #buffers
+        mf = cl.mem_flags
+        self.qx_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
+        self.qy_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
+        self.res_b = cl.Buffer(self.ctx, mf.WRITE_ONLY, qx.nbytes)
+        self.res = np.empty_like(self.qx)
+
+    def eval(self, pars):
+
+        bi_thick = GaussianDispersion(int(pars['bi_thick_pd_n']), pars['bi_thick_pd'], pars['bi_thick_pd_nsigma'])
+        bi_thick.value, bi_thick.weight = bi_thick.get_weights(pars['bi_thick'], 0, 1000, True)
+
+        sum, norm = 0.0, 0.0
+        sub = pars['sld_bi'] - pars['sld_sol']
+
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
+        for i in xrange(len(bi_thick.weight)):
+            self.prg.LamellarKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(bi_thick.value[i]),
+                                    real(pars['scale']), real(sub), np.uint32(self.qx.size))
+            cl.enqueue_copy(self.queue, self.res, self.res_b)
+
+            sum += bi_thick.weight[i]*self.res
+            norm += bi_thick.weight[i]
+
+        return sum/norm + pars['background']
+
+def demo():
+
+    from time import time
+    import matplotlib.pyplot as plt
+
+    #create qx and qy evenly spaces
+    qx = np.linspace(-.01, .01, 128)
+    qy = np.linspace(-.01, .01, 128)
+    qx, qy = np.meshgrid(qx, qy)
+
+    #saved shape of qx
+    r_shape = qx.shape
+
+    #reshape for calculation; resize as float32
+    qx = qx.flatten()
+    qy = qy.flatten()
+
+    pars = LamellarParameters(scale=1, bi_thick=100, sld_bi=.291e-6, sld_sol=5.77e-6, background=0)
+    t = time()
+    result = GpuLamellar(qx, qy)
+    result.x = result.lamellar_fit(pars, b_n=35, b_w=.1, sigma=3)
+    result.x = np.reshape(result.x, r_shape)
+    tt = time()
+
+    print("Time taken: %f" % (tt - t))
+
+    f = open("r.txt", "w")
+    for x in xrange(len(r_shape)):
+        f.write(str(result.x[x]))
+        f.write("\n")
+
+    plt.pcolormesh(np.log10(result.x))
+    plt.show()
+
+
+if __name__ == "__main__":
+    demo()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: code_triaxialellipse.py
===================================================================
--- code_triaxialellipse.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
+++ code_triaxialellipse.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -0,0 +1,128 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+from math import asin
+import pyopencl as cl
+from weights import GaussianDispersion
+from sasmodel import card
+
+def set_precision(src, qx, qy, dtype):
+    qx = np.ascontiguousarray(qx, dtype=dtype)
+    qy = np.ascontiguousarray(qy, dtype=dtype)
+    if np.dtype(dtype) == np.dtype('float32'):
+        header = """\
+#define real float
+"""
+    else:
+        header = """\
+#pragma OPENCL EXTENSION cl_khr_fp64: enable
+#define real double
+"""
+    return header+src, qx, qy
+
+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}
+
+    PD_PARS = ['axisA', 'axisB', 'axisC', 'theta', 'phi', 'psi']
+
+    def __init__(self, qx, qy, dtype='float32'):
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('Kernel-TriaxialEllipse.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+        #buffers
+        mf = cl.mem_flags
+        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.WRITE_ONLY, qx.nbytes)
+        self.res = np.empty_like(self.qx)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        axisA, axisB, axisC, theta, phi, 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, 1000, True)
+        axisB.value, axisB.weight = axisB.get_weights(pars['axisB'], 0, 1000, True)
+        axisC.value, axisC.weight = axisC.get_weights(pars['axisC'], 0, 1000, 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)
+
+        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
+        size = len(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)):
+                                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))
+                                cl.enqueue_copy(queue, self.res, self.res_b)
+                                sum += self.res
+
+                                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]
+
+        if size > 1:
+            norm /= asin(1.0)
+
+        if vol != 0.0 and norm_vol != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm + pars['background']
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: reshellcylcode.py
===================================================================
--- coreshellcylcode.py	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ 	(revision )
@@ -1,86 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import math
-import pyopencl as cl
-from weights import GaussianDispersion
-from sasmodel import card
-
-def set_precision(src, qx, qy, dtype):
-    qx = np.ascontiguousarray(qx, dtype=dtype)
-    qy = np.ascontiguousarray(qy, dtype=dtype)
-    if np.dtype(dtype) == np.dtype('float32'):
-        header = """\
-#define real float
-"""
-    else:
-        header = """\
-#pragma OPENCL EXTENSION cl_khr_fp64: enable
-#define real double
-"""
-    return header+src, qx, qy
-
-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'):
-        #create context, queue, and build program
-        ctx,_queue = card()
-        src, qx, qy = set_precision(open('NR_BessJ1.cpp').read()+"\n"+open('Kernel-CoreShellCylinder.cpp').read(), qx, qy, dtype=dtype)
-        self.prg = cl.Program(ctx, src).build()
-        self.qx, self.qy = qx, qy
-
-
-        #buffers
-        mf = cl.mem_flags
-        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.WRITE_ONLY, qx.nbytes)
-        self.res = np.empty_like(qx)
-
-    def eval(self, pars):
-
-        _ctx,queue = card()
-        radius, length, thickness, axis_phi, axis_theta = [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
-                                     for base in GpuCoreShellCylinder.PD_PARS]
-
-        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 1000, True)
-        length.value, length.weight = length.get_weights(pars['length'], 0, 1000, True)
-        thickness.value, thickness.weight = thickness.get_weights(pars['thickness'], 0, 1000, True)
-        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False)
-        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False)
-
-        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
-        size = len(axis_theta.weight)
-
-        for i in xrange(len(radius.weight)):
-            for j in xrange(len(length.weight)):
-                for k in xrange(len(axis_theta.weight)):
-                    for l in xrange(len(axis_phi.weight)):
-                        for f in xrange(len(thickness.weight)):
-
-                            self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                    np.float32(axis_theta.value[k]), np.float32(axis_phi.value[l]), np.float32(thickness.value[f]),
-                                    np.float32(length.value[j]), np.float32(radius.value[i]), np.float32(pars['scale']),
-                                    np.float32(radius.weight[i]), np.float32(length.weight[j]), np.float32(thickness.weight[f]),
-                                    np.float32(axis_theta.weight[k]), np.float32(axis_phi.weight[l]), np.float32(pars['core_sld']),
-                                    np.float32(pars['shell_sld']), np.float32(pars['solvent_sld']),np.uint32(size),
-                                    np.uint32(self.qx.size))
-                            cl.enqueue_copy(queue, self.res, self.res_b)
-
-                            sum += self.res
-                            vol += radius.weight[i]*length.weight[j]*thickness.weight[f]*pow(radius.value[i]+thickness.value[f],2)\
-                                   *(length.value[j]+2.0*thickness.value[f])
-                            norm_vol += radius.weight[i]*length.weight[j]*thickness.weight[k]
-                            norm += radius.weight[i]*length.weight[j]*thickness.weight[f]*axis_theta.weight[k]\
-                                    *axis_phi.weight[l]
-
-        if size>1:
-            norm /= math.asin(1.0)
-        if vol != 0.0 and norm_vol != 0.0:
-            sum *= norm_vol/vol
-
-        return sum/norm + pars['background']
Index: lcode.py
===================================================================
--- cylcode.py	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ 	(revision )
@@ -1,148 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import math
-import pyopencl as cl
-from weights import GaussianDispersion
-from sasmodel import card
-
-def set_precision(src, qx, qy, dtype):
-    qx = np.ascontiguousarray(qx, dtype=dtype)
-    qy = np.ascontiguousarray(qy, dtype=dtype)
-    if np.dtype(dtype) == np.dtype('float32'):
-        header = """\
-#define real float
-"""
-    else:
-        header = """\
-#pragma OPENCL EXTENSION cl_khr_fp64: enable
-#define real double
-"""
-    return header+src, qx, qy
-
-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'):
-
-        #create context, queue, and build program
-        ctx,_queue = card()
-        src, qx, qy = set_precision(open('NR_BessJ1.cpp').read()+"\n"+open('Kernel-Cylinder.cpp').read(), qx, qy, dtype=dtype)
-        self.prg = cl.Program(ctx, src).build()
-        self.qx, self.qy = qx, qy
-
-        #buffers
-        mf = cl.mem_flags
-        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.WRITE_ONLY, qx.nbytes)
-        self.res = np.empty_like(self.qx)
-
-    def eval(self, pars):
-
-        _ctx,queue = card()
-        radius, length, cyl_theta, cyl_phi = \
-            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
-             for base in GpuCylinder.PD_PARS]
-
-        #Get the weights for each
-        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 1000, True)
-        length.value, length.weight = length.get_weights(pars['length'], 0, 1000, True)
-        cyl_theta.value, cyl_theta.weight = cyl_theta.get_weights(pars['cyl_theta'], -90, 180, False)
-        cyl_phi.value, cyl_phi.weight = cyl_phi.get_weights(pars['cyl_phi'], -90, 180, False)
-
-        #Perform the computation, with all weight points
-        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
-        size = len(cyl_theta.weight)
-        sub = pars['sldCyl'] - pars['sldSolv']
-
-        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-        #Loop over radius, length, theta, phi weight points
-        for i in xrange(len(radius.weight)):
-            for j in xrange(len(length.weight)):
-                for k in xrange(len(cyl_theta.weight)):
-                    for l in xrange(len(cyl_phi.weight)):
-
-
-                        self.prg.CylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(sub),
-                                           real(radius.value[i]), real(length.value[j]), real(pars['scale']),
-                                           real(radius.weight[i]), real(length.weight[j]), real(cyl_theta.weight[k]),
-                                           real(cyl_phi.weight[l]), real(cyl_theta.value[k]), real(cyl_phi.value[l]),
-                                           np.uint32(self.qx.size), np.uint32(size))
-                        cl.enqueue_copy(queue, self.res, self.res_b)
-                        sum += self.res
-                        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:
-            norm /= math.asin(1.0)
-        if vol != 0.0 and norm_vol != 0.0:
-            sum *= norm_vol/vol
-
-        return sum/norm+pars['background']
-
-def demo():
-    from time import time
-    import matplotlib.pyplot as plt
-
-    #create qx and qy evenly spaces
-    qx = np.linspace(-.02, .02, 128)
-    qy = np.linspace(-.02, .02, 128)
-    qx, qy = np.meshgrid(qx, qy)
-
-    #saved shape of qx
-    r_shape = qx.shape
-
-    #reshape for calculation; resize as float32
-    qx = qx.flatten()
-    qy = qy.flatten()
-
-    pars = CylinderParameters(scale=1, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=0,
-                              cyl_theta=0, cyl_phi=0)
-    t = time()
-    result = GpuCylinder(qx, qy)
-    result.x = result.cylinder_fit(pars, r_n=10, t_n=10, l_n=10, p_n=10, r_w=.1, t_w=.1, l_w=.1, p_w=.1, sigma=3)
-    result.x = np.reshape(result.x, r_shape)
-    tt = time()
-
-    print("Time taken: %f" % (tt - t))
-
-    plt.pcolormesh(result.x)
-    plt.show()
-
-
-if __name__=="__main__":
-    demo()
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: lipsecode.py
===================================================================
--- ellipsecode.py	(revision 5378e40a7fb092abdefdce6d59001232f68a74e2)
+++ 	(revision )
@@ -1,113 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import math
-import pyopencl as cl
-from weights import GaussianDispersion
-
-
-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):
-
-        self.qx = np.asarray(qx, np.float32)
-        self.qy = np.asarray(qy, np.float32)
-        #create context, queue, and build program
-        self.ctx = cl.create_some_context()
-        self.queue = cl.CommandQueue(self.ctx)
-        self.prg = cl.Program(self.ctx, open('Kernel-Ellipse.cpp').read()).build()
-
-        #buffers
-        mf = cl.mem_flags
-        self.qx_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
-        self.qy_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
-        self.res_b = cl.Buffer(self.ctx, mf.WRITE_ONLY, qx.nbytes)
-        self.res = np.empty_like(self.qx)
-
-    def eval(self, pars):
-    #b_n = radius_b # want, a_n = radius_a # want, etc
-        radius_a, radius_b, axis_theta, axis_phi = \
-            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
-             for base in GpuEllipse.PD_PARS]
-
-        radius_a.value, radius_a.weight = radius_a.get_weights(pars['radius_a'], 0, 1000, True)
-        radius_b.value, radius_b.weight = radius_b.get_weights(pars['radius_b'], 0, 1000, 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)
-
-        #Perform the computation, with all weight points
-        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
-        size = len(axis_theta.weight)
-        sub =  pars['sldEll'] - pars['sldSolv']
-
-        #Loop over radius weight points
-        for i in xrange(len(radius_a.weight)):
-            #Loop over length weight points
-            for j in xrange(len(radius_b.weight)):
-                #Average over theta distribution
-                for k in xrange(len(axis_theta.weight)):
-                    #Average over phi distribution
-                    for l in xrange(len(axis_phi.weight)):
-                        #call the kernel
-                        self.prg.EllipsoidKernel(self.queue, self.qx.shape, None, np.float32(radius_a.weight[i]),
-                                        np.float32(radius_b.weight[j]), np.float32(axis_theta.weight[k]),
-                                        np.float32(axis_phi.weight[l]), np.float32(pars['scale']), np.float32(radius_a.value[i]),
-                                        np.float32(radius_b.value[j]), np.float32(sub),np.float32(pars['background']),
-                                        np.float32(axis_theta.value[k]), np.float32(axis_phi.value[l]), self.qx_b, self.qy_b,
-                                        self.res_b, np.uint32(self.qx.size), np.uint32(len(axis_theta.weight)))
-                        #copy result back from buffer
-                        cl.enqueue_copy(self.queue, self.res, self.res_b)
-                        sum += self.res
-                        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
-        # integration (see documentation).
-        if size > 1:
-            norm /= math.asin(1.0)
-        if vol.any() != 0.0 and norm_vol.any() != 0.0:
-            sum *= norm_vol/vol
-
-        return sum/norm+pars['background']
-
-
-def demo():
-    from time import time
-    import matplotlib.pyplot as plt
-
-    #create qx and qy evenly spaces
-    qx = np.linspace(-.02, .02, 128)
-    qy = np.linspace(-.02, .02, 128)
-    qx, qy = np.meshgrid(qx, qy)
-
-    #saved shape of qx
-    r_shape = qx.shape
-    #reshape for calculation; resize as float32
-    qx = qx.flatten()
-    qy = qy.flatten()
-
-    #int main
-    pars = EllipsoidParameters(.027, 60, 180, .297e-6, 5.773e-06, 4.9, 0, 90)
-
-    t = time()
-    result = GpuEllipse(qx, qy)
-    result.x = result.ellipsoid_fit(qx, qy, pars, b_n=35, t_n=35, a_n=1, p_n=1, sigma=3, b_w=.1, t_w=.1, a_w=.1, p_w=.1)
-    result.x = np.reshape(result.x, r_shape)
-    tt = time()
-    print("Time taken: %f" % (tt - t))
-
-    plt.pcolormesh(result.x)
-    plt.show()
-
-
-if __name__ == "__main__":
-    demo()
-
-
-
-
Index: fit.py
===================================================================
--- fit.py	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
+++ fit.py	(revision 2de9a5e5fa221a201e0d071bc49f8bf651f6f679)
@@ -3,8 +3,10 @@
 
 from bumps.names import *
-from cylcode import GpuCylinder
-from lamellarcode import GpuLamellar
-from ellipsecode import GpuEllipse
-from coreshellcylcode import GpuCoreShellCylinder
+from code_cylinder import GpuCylinder
+from code_lamellar import GpuLamellar
+from code_ellipse import GpuEllipse
+from code_coreshellcyl import GpuCoreShellCylinder
+from code_capcyl import GpuCapCylinder
+from code_triaxialellipse import GpuTriEllipse
 from sasmodel import SasModel, load_data, set_beam_stop
 
@@ -25,20 +27,28 @@
                  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, dtype='float')
+                 axis_phi_pd_n=6, axis_phi_pd_nsigma=3, dtype='float32')
 
 model = SasModel(data, GpuLamellar, scale=1, bi_thick=100, sld_bi=.291e-6, sld_sol=5.77e-6, background=0,
-                 bi_thick_pd=0.1, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float')
+                 bi_thick_pd=0.1, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float32')
 
-"""
+
 model = SasModel(data, GpuCoreShellCylinder, scale=1, radius=64.1, thickness=1, length=266.96, core_sld=1e-6, shell_sld=4e-6,
                  solvent_sld=1e-6, background=0, axis_theta=0, 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=2, thickness_pd_nsigma=3,
                  axis_theta_pd=0.1, axis_theta_pd_n=2, axis_theta_pd_nsigma=3, axis_phi_pd=0.1, axis_phi_pd_n=2,
-                 axis_phi_pd_nsigma=3, dtype='float')
+                 axis_phi_pd_nsigma=3, dtype='float32')
 
-model = SasModel(data, GpuCapCylinder, scale, rad_cyl, rad_cap, length, sld_capcyl, sld_solv, background, theta, phi,
-                 rad_cyl_pd, rad_cyl_pd_n, rad_cyl_nsigma, rad_cap_pd, rad_cap_pd_n, rad_cap_pd_nsigma, length_)
+model = SasModel(data, GpuCapCylinder, scale=1, rad_cyl=20, rad_cap=40, length=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_nsigma=3, rad_cap_pd=.1, rad_cap_pd_n=1,
+                 rad_cap_pd_nsigma=3, length_pd=.1, length_pd_n=10, length_pd_nsigma=3, theta_pd=.1, theta_pd_n=4,
+                 theta_pd_nsigma=3, phi_pd=.1, phi_pd_n=4, phi_pd_nsigma=3, dtype='float32')
+"""
+model = SasModel(data, GpuTriEllipse, scale=1, axisA=35, axisB=100, axisC=400, sldEll=1e-6, sldSolv=6.3e-6,
+                 background=0, theta=57, phi=57, psi=0, theta_pd=.1, theta_pd_n=4,
+                 theta_pd_nsigma=3, phi_pd=.1, phi_pd_n=4, phi_pd_nsigma=3, psi_pd=.1, psi_pd_n=4, psi_pd_nsigma=3,
+                 axisA_pd=.1, axisA_pd_n=4, axisA_pd_nsigma=3, axisB_pd=.1, axisB_pd_n=4, axisB_pd_nsigma=3,
+                 axisC_pd=.1, axisC_pd_n=4, axisC_pd_nsigma=3, dtype='float')
 
-model.scale.range(0,10)
+model.scale.range(0,1)
 
 problem = FitProblem(model)
Index: mellarcode.py
===================================================================
--- lamellarcode.py	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ 	(revision )
@@ -1,129 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import pyopencl as cl
-from weights import GaussianDispersion
-
-def set_precision(src, qx, qy, dtype):
-    qx = np.ascontiguousarray(qx, dtype=dtype)
-    qy = np.ascontiguousarray(qy, dtype=dtype)
-    if dtype == 'double':
-        header = """\
-#pragma OPENCL EXTENSION cl_khr_fp64: enable
-#define real double
-"""
-        return header+src,qx,qy
-    else:
-        return src,qx,qy
-
-
-class GpuLamellar(object):
-    PARS = {
-        'scale':1, 'bi_thick':1, 'sld_bi':1e-6, 'sld_sol':0, 'background':0,
-    }
-
-    def __init__(self, qx, qy, dtype='float'):
-
-        #create context, queue, and build program
-        self.ctx = cl.create_some_context()
-        self.queue = cl.CommandQueue(self.ctx)
-        src,qx,qy = set_precision(open('Kernel-Lamellar.cpp').read(), qx, qy, dtype=dtype)
-        self.prg = cl.Program(self.ctx, src).build()
-        self.qx, self.qy = qx, qy
-
-        #buffers
-        mf = cl.mem_flags
-        self.qx_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
-        self.qy_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
-        self.res_b = cl.Buffer(self.ctx, mf.WRITE_ONLY, qx.nbytes)
-        self.res = np.empty_like(self.qx)
-
-    def eval(self, pars):
-
-        bi_thick = GaussianDispersion(int(pars['bi_thick_pd_n']), pars['bi_thick_pd'], pars['bi_thick_pd_nsigma'])
-        bi_thick.value, bi_thick.weight = bi_thick.get_weights(pars['bi_thick'], 0, 1000, True)
-
-        sum, norm = 0.0, 0.0
-        sub = pars['sld_bi'] - pars['sld_sol']
-
-        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-        for i in xrange(len(bi_thick.weight)):
-            self.prg.LamellarKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(bi_thick.value[i]),
-                                    real(pars['scale']), real(sub), real(pars['background']), np.uint32(self.qx.size))
-            cl.enqueue_copy(self.queue, self.res, self.res_b)
-
-            sum += bi_thick.weight[i]*self.res
-            norm += bi_thick.weight[i]
-
-        return sum/norm + pars['background']
-
-def demo():
-
-    from time import time
-    import matplotlib.pyplot as plt
-
-    #create qx and qy evenly spaces
-    qx = np.linspace(-.01, .01, 128)
-    qy = np.linspace(-.01, .01, 128)
-    qx, qy = np.meshgrid(qx, qy)
-
-    #saved shape of qx
-    r_shape = qx.shape
-
-    #reshape for calculation; resize as float32
-    qx = qx.flatten()
-    qy = qy.flatten()
-
-    pars = LamellarParameters(scale=1, bi_thick=100, sld_bi=.291e-6, sld_sol=5.77e-6, background=0)
-    t = time()
-    result = GpuLamellar(qx, qy)
-    result.x = result.lamellar_fit(pars, b_n=35, b_w=.1, sigma=3)
-    result.x = np.reshape(result.x, r_shape)
-    tt = time()
-
-    print("Time taken: %f" % (tt - t))
-
-    f = open("r.txt", "w")
-    for x in xrange(len(r_shape)):
-        f.write(str(result.x[x]))
-        f.write("\n")
-
-    plt.pcolormesh(np.log10(result.x))
-    plt.show()
-
-
-if __name__ == "__main__":
-    demo()
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
