Index: Capcyl_Kfun.cpp
===================================================================
--- Capcyl_Kfun.cpp	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
+++ Capcyl_Kfun.cpp	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
@@ -0,0 +1,171 @@
+real ConvLens_kernel(real len, real rad, real endRad, real x, real tt, real theta)
+{
+    real be = 0;
+	real hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));
+	real arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
+	real arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
+
+	if(arg2 == 0) {be = 0.5;}
+	else {
+		be = NR_BessJ1(arg2)/arg2;
+	}
+	real val = cos(arg1)*(1.0-tt*tt)*be;
+
+	return(val);
+}
+const real Gauss76Z[76]={
+	 .999505948362153*(-1.0),		//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),		//10
+	 .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),		//20
+	 .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),		//30
+	 .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),			//40
+	-.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),		//50
+	-.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),	//60
+	-.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),		//70
+	-.981013938975656*(-1.0),
+	-.988144453359837*(-1.0),
+	-.993608772723527*(-1.0),
+	-.997397786355355*(-1.0),
+	-.999505948362153*(-1.0)		//75
+};
+const real Gauss76Wt[76]={
+	.00126779163408536,		//0
+	.00294910295364247,
+	.00462793522803742,
+	.00629918049732845,
+	.00795984747723973,
+	.00960710541471375,
+	.0112381685696677,
+	.0128502838475101,
+	.0144407317482767,
+	.0160068299122486,
+	.0175459372914742,		//10
+	.0190554584671906,
+	.020532847967908,
+	.0219756145344162,
+	.0233813253070112,
+	.0247476099206597,
+	.026072164497986,
+	.0273527555318275,
+	.028587223650054,
+	.029773487255905,
+	.0309095460374916,		//20
+	.0319934843404216,
+	.0330234743977917,
+	.0339977794120564,
+	.0349147564835508,
+	.0357728593807139,
+	.0365706411473296,
+	.0373067565423816,
+	.0379799643084053,
+	.0385891292645067,
+	.0391332242205184,		//30
+	.0396113317090621,
+	.0400226455325968,
+	.040366472122844,
+	.0406422317102947,
+	.0408494593018285,
+	.040987805464794,
+	.0410570369162294,
+	.0410570369162294,
+	.040987805464794,
+	.0408494593018285,		//40
+	.0406422317102947,
+	.040366472122844,
+	.0400226455325968,
+	.0396113317090621,
+	.0391332242205184,
+	.0385891292645067,
+	.0379799643084053,
+	.0373067565423816,
+	.0365706411473296,
+	.0357728593807139,		//50
+	.0349147564835508,
+	.0339977794120564,
+	.0330234743977917,
+	.0319934843404216,
+	.0309095460374916,
+	.029773487255905,
+	.028587223650054,
+	.0273527555318275,
+	.026072164497986,
+	.0247476099206597,		//60
+	.0233813253070112,
+	.0219756145344162,
+	.020532847967908,
+	.0190554584671906,
+	.0175459372914742,
+	.0160068299122486,
+	.0144407317482767,
+	.0128502838475101,
+	.0112381685696677,
+	.00960710541471375,		//70
+	.00795984747723973,
+	.00629918049732845,
+	.00462793522803742,
+	.00294910295364247,
+	.00126779163408536		//75 (indexed from 0)
+};
Index: Kernel-CapCyl.cpp
===================================================================
--- Kernel-CapCyl.cpp	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ Kernel-CapCyl.cpp	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
@@ -22,8 +22,8 @@
         float vaj = -1.0*hDist/rad_cap;
 
-        for(j=0;j<76;j++) //the value 76 was pre-set in the kernel...
+        for(j=0;j<76;j++) //the 76 corresponds to the Gauss constants
         {
-            zij = ( Gauss76Z[j]*(1.0-vaj) + vaj + 1.0 )/2.0;    //the "t" dummy
-            yyy = Gauss76Wt[j]*ConvLens_kernel(dp,q,zij,alpha);    //uses the same Kernel as the Dumbbell, here L>0
+            zij = (Gauss76Z[j]*(1.0-vaj)+vaj+1.0)/2.0;
+            yyy = Gauss76Wt[j]*ConvLens_kernel(length,rad_cyl,rad_cap,q,zij,alpha);
             summj += yyy;
         }
@@ -35,23 +35,8 @@
         if(arg2 == 0) {be = 0.5;}
         else {
-            float ax=fabs(arg2);
-            if ((ax < 8.0) {
-                y=arg2*arg2;
-                ans1=arg2*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
-                ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))));
-                ans=ans1/ans2;
-            } 
-            else {
-                y=64.0/(ax*ax);
-                xx=ax-2.356194491;
-                ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6))));
-                ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)));
-                ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-(8.0/ax)*sin(xx)*ans2);
-                if (arg2 < 0.0) {ans *= -1;}
-	        }
-	        be = ans/arg2
+            be = NR_BessJ1(arg2)/arg2;
         }
 
-        if(arg1 == 0.0) {   //limiting value of sinc(0) is 1; sinc is not defined in math.h
+        if(arg1 == 0.0) {
             yyy += pi*rad_cyl*rad_cyl*length*2.0*be;
         } 
Index: Kernel-CoreShellCylinder.cpp
===================================================================
--- Kernel-CoreShellCylinder.cpp	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ Kernel-CoreShellCylinder.cpp	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
@@ -22,5 +22,4 @@
 
         float si1=0; float si2=0; float be1=0; float be2=0;
-        float ax=0; float z=0; float xx=0; float y=0; float ans=0; float ans1=0; float ans2=0;
 
         float dr1 = core_sld-shell_sld;
@@ -37,26 +36,9 @@
         if (besarg1 == 0.0){be1 = 0.5;}
         else{
-            bel = NR_BessJ1(besarg1)/besarg1
+            be1 = NR_BessJ1(besarg1)/besarg1;
         }
         if (besarg2 == 0.0){be2 = 0.5;}
         else{
-            if((ax=fabs(besarg2)) < 8.0)
-            {
-                y=besarg2*besarg2;
-                ans1=besarg2*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
-                ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))));
-                ans=ans1/ans2;
-            }
-            else
-            {
-                z=8.0/ax;
-                y=z*z;
-                xx=ax-2.356194491;
-                ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6))));
-                ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)));
-                ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
-                if (besarg2 < 0.0) {ans *= -1;}
-            }
-            be2 = ans/besarg2;
+            be2 = NR_BessJ1(besarg2)/besarg2;
         }
         if (sinarg1 == 0.0){
Index: NR_BessJ1.cpp
===================================================================
--- NR_BessJ1.cpp	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ NR_BessJ1.cpp	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
Index: capcylcope.py
===================================================================
--- capcylcope.py	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ capcylcope.py	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
@@ -6,4 +6,19 @@
 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):
@@ -13,25 +28,25 @@
     PD_PARS = ['rad_cyl', 'length', 'rad_cap', 'theta', 'phi']
 
-    def __init__(self, qx, qy):
+    def __init__(self, qx, qy, dtype='float32'):
 
-        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)
+        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(self.ctx, open('Kernel-CapCyl.cpp').read()).build()
+        self.prg = cl.Program(ctx, open('Kernel-CapCyl.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.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(self.ctx, mf.WRITE_ONLY, self.vol_i.nbytes)
+        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'])
@@ -54,5 +69,5 @@
                         for l in xrange(len(phi.weight)):
 
-                            self.prg.CapCylinderKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
+                            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']),
@@ -60,6 +75,6 @@
                                         np.float32(rad_cyl.weight[i]), np.float32(length.weight[j]), np.uint32(self.qx.size), np.uint32(size))
 
-                            cl.enqueue_copy(self.queue, self.res, self.res_b)
-                            cl.enqueue_copy(self.queue, self.vol_i, self.vol_b)
+                            cl.enqueue_copy(queue, self.res, self.res_b)
+                            cl.enqueue_copy(queue, self.vol_i, self.vol_b)
 
                             sum += self.res
Index: fit.py
===================================================================
--- fit.py	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ fit.py	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
@@ -12,8 +12,5 @@
 data = load_data('JUN03289.DAT')
 set_beam_stop(data, 0.004)
-
-
-
-
+"""
 model = SasModel(data, GpuCylinder, scale=1, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=0,
                               cyl_theta=0, cyl_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,length_pd=0.1,
@@ -24,5 +21,5 @@
 model.cyl_theta.range(0,90)
 model.cyl_phi.range(0,90)
-"""
+
 model = SasModel(data, GpuEllipse, 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,
@@ -33,11 +30,13 @@
                  bi_thick_pd=0.1, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float')
 
-
-model = SasModel(data, GpuCoreShellCylinder, scale=1, radius=64.1, thickness=1, length=266.96, core_sld=1e-6, shell_sld=1e-6,
-                 solvent_sld=4e-6, background=0, axis_theta=0, axis_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
+"""
+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')
-"""
+
+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.scale.range(0,10)
Index: fit2.py
===================================================================
--- fit2.py	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ fit2.py	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
@@ -34,10 +34,8 @@
 radial.cyl_phi.range(0,90)
 radial.scale.range(0,10)
+trans.cyl_theta = radial.cyl_theta + 90.
 
 
-fv = FreeVariables(names=['radial','trans'],
-                   cyl_theta=radial.cyl_theta,
-                   cyl_phi = radial.cyl_phi)
-problem = FitProblem([radial,trans], freevars=fv)
+problem = FitProblem([radial,trans])
 
 
Index: sasmodel.py
===================================================================
--- sasmodel.py	(revision 8a20be51c8c1b5b0757f10d24eb2bba0e3ddff6c)
+++ sasmodel.py	(revision 8faffcd2c614126a3e6035d58f1f258f3271103e)
@@ -62,4 +62,5 @@
 class SasModel(object):
     def __init__(self, data, model, dtype='float32', **kw):
+        self.__dict__['_parameters'] = {}
         self.index = data.mask==0
         self.iq = data.data[self.index]
@@ -79,5 +80,5 @@
         pars.update((base+'_pd_nsigma', 3) for base in model.PD_PARS)
         pars.update(kw)
-        self._parameters = dict((k, Parameter(v, name=k)) for k, v in pars.items())
+        self._parameters = dict((k, Parameter.default(v, name=k)) for k, v in pars.items())
 
     def numpoints(self):
@@ -90,7 +91,12 @@
         return self._parameters[par]
 
+    def __setattr__(self, par, val):
+        if par in self._parameters:
+            self._parameters[par] = val
+        else:
+            self.__dict__[par] = val
+
     def theory(self):
         pars = dict((k,v.value) for k,v in self._parameters.items())
-        print pars
         result = self.gpu.eval(pars)
         return result
