Index: rnel/Capcyl_Kfun.cpp
===================================================================
--- Kernel/Capcyl_Kfun.cpp	(revision 099e05357033bb53919d7a0213dc1bd7a0b54ebb)
+++ 	(revision )
@@ -1,63 +1,0 @@
-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);
-}
-real Gauss76Z(int i)
-{
-real array[76] = {.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) };
-
-return array[i];
-}
-real Gauss76Wt(int i)
-{
-real array[76] = {.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 };
-
-return array[i];
-}
Index: rnel/Kernel-CapCyl.cpp
===================================================================
--- Kernel/Kernel-CapCyl.cpp	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ 	(revision )
@@ -1,49 +1,0 @@
-__kernel void CapCylinderKernel(__global const real *qx, __global const real *qy, __global real *_ptvalue, const real vol_i,
-const real hDist, const real rad_cyl, const real rad_cap, const real length, const real thet, const real ph, const real sub,
-const real scale, const real phi_weight, const real theta_float, const real rad_cap_weight, const real rad_cyl_weight,
-const real length_weight, const real theta_weight, const int total, const int size)
-//ph is phi, sub is sldc-slds, thet is theta
-{
-    int i = get_global_id(0);
-    if(i < total)
-    {
-        real q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
-        real pi = 4.0*atan(1.0);
-        real theta = thet*pi/180.0;
-        real phi = ph*pi/180.0;
-        real cyl_x = cos(theta)*cos(phi);
-        real cyl_y = sin(theta);
-        real cos_val = cyl_x*qx[i]/q + cyl_y*qy[i]/q;
-        real alpha = acos(cos_val);
-        real ans1=0; real ans2=0; real y=0; real xx=0; real ans=0; real zij=0; real be=0; real summj=0;
-        real vaj = -1.0*hDist/rad_cap;
-
-        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;
-            summj += Gauss76Wt(j)*ConvLens_kernel(length,rad_cyl,rad_cap,q,zij,alpha);
-        }
-        real yyy = (1.0-vaj)/2.0*summj*4.0*pi*rad_cap*rad_cap*rad_cap;
-        real arg1 = q*length/2.0*cos(alpha);
-        real arg2 = q*rad_cyl*sin(alpha);
-
-        if(arg2 == 0) {be = 0.5;}
-        else {
-            be = NR_BessJ1(arg2)/arg2;
-        }
-
-        if(arg1 == 0.0) {
-            yyy += pi*rad_cyl*rad_cyl*length*2.0*be;
-        } 
-        else {
-            yyy += pi*rad_cyl*rad_cyl*length*sin(arg1)/arg1*2.0*be;
-        }
-        real 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);
-
-        _ptvalue[i] += rad_cyl_weight*length_weight*rad_cap_weight*theta_weight*phi_weight*vol_i*answer;
-     //   if (size>1) {
-   //         _ptvalue[i] *= fabs(cos(thet*pi/180.0));
-     //   }
-    }
-}
Index: rnel/Kernel-CoreShellCylinder.cpp
===================================================================
--- Kernel/Kernel-CoreShellCylinder.cpp	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ 	(revision )
@@ -1,55 +1,0 @@
-__kernel void CoreShellCylinderKernel(__global const real *qx, __global const real *qy, __global real *_ptvalue,
-const real axis_theta, const real axis_phi, const real thickness, const real length, const real radius,
-const real scale, const real radius_weight, const real length_weight, const real thickness_weight,
-const real theta_weight, const real phi_weight, const real core_sld, const real shell_sld, const real solvent_sld,
-const int size, const int total)
-{
-    int i = get_global_id(0);
-    if(i < total)
-    {
-        real q = sqrt(qx[i]*qx[i]+qy[i]*qy[i]);
-        real pi = 4.0*atan(1.0);
-        real theta = axis_theta*pi/180.0;
-        real alpha = acos(cos(theta)*cos(axis_phi*pi/180.0)*qx[i]/q + sin(theta)*qy[i]/q);
-
-        if (alpha == 0.0){
-        alpha = 1.0e-26;
-        }
-
-        real si1=0; real si2=0; real be1=0; real be2=0;
-
-        real vol2 = pi*(radius+thickness)*(radius+thickness)*(length+2.0*thickness);
-
-        real besarg1 = q*radius*sin(alpha);
-        real besarg2 = q*(radius+thickness)*sin(alpha);
-        real sinarg1 = q*length/2.0*cos(alpha);
-        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;}
-
-        real tt = 2.0*vol2*(shell_sld-solvent_sld)*si2*be2+2.0*(pi*radius*radius*(length))*(core_sld-shell_sld)*si1*be1;
-
-        real answer = (tt*tt)*sin(alpha)/fabs(sin(alpha));
-        answer *= answer/(pi*(radius+thickness)*(radius+thickness)*(length+2.0*thickness))*1.0e8*scale;
-
-        _ptvalue[i] += radius_weight*length_weight*thickness_weight*theta_weight*phi_weight*answer*pow(radius+thickness,2)*(length+2.0*thickness);
-     //   if (size>1) {
-       // _ptvalue[i] *= fabs(cos(axis_theta*pi/180.0));
-        //}
-
-    }
-}
-
-
-
Index: rnel/Kernel-CoreShellCylinder_f.cpp
===================================================================
--- Kernel/Kernel-CoreShellCylinder_f.cpp	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
+++ 	(revision )
@@ -1,128 +1,0 @@
-#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: rnel/Kernel-Cylinder.cpp
===================================================================
--- Kernel/Kernel-Cylinder.cpp	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ 	(revision )
@@ -1,79 +1,0 @@
-__kernel void CylinderKernel(__global const real *qx, global const real *qy, __global real *_ptvalue, const real sub,
-const real rr, const real h, const real scale, const real radius_weight, const real length_weight,
-const real theta_weight, const real phi_weight, const real cyl_theta,
-const real cyl_phi, const int count, const int size)
-{
-	// qq is the q-value for the calculation (1/A)
-	// rr is the radius of the cylinder (A)
-	// h is the -LENGTH of the cylinder = L (A)
-    int i = get_global_id(0);
-
-    if(i < count)
-    {
-        real be=0.0; real si=0.0;
-        real qq = sqrt(qx[i]*qx[i]+qy[i]*qy[i]);
-
-        real pi = 4.0*atan(1.0);
-        real theta = cyl_theta*pi/180.0;
-        real cos_val = cos(theta)*cos(cyl_phi*pi/180.0)*(qx[i]/qq) + sin(theta)*(qy[i]/qq);
-
-        real alpha = acos(cos_val);
-        if(alpha == 0.0){
-            alpha = 1.0e-26;
-        }
-        real besarg = qq*rr*sin(alpha);
-        real siarg = qq*h/2*cos(alpha);
-
-        real bj = NR_BessJ1(besarg);
-        real d1 = qq*rr*sin(alpha);
-
-        if (besarg == 0.0){
-            be = sin(alpha);
-        }
-        else{
-            be = bj*bj*4.0*sin(alpha)/(d1*d1);
-        }
-        if(siarg == 0.0){
-            si = 1.0;
-        }else{
-            si = sin(siarg)*sin(siarg)/(siarg*siarg);
-        }
-        real form = be*si/sin(alpha);
-        real answer = sub*sub*form*acos(-1.0)*rr*rr*h*1.0e8*scale;
-
-        _ptvalue[i] += radius_weight*length_weight*theta_weight*phi_weight*answer*pow(rr,2)*h;
-
-       // if (size>1) {
-         //   _ptvalue[i] *= fabs(cos(cyl_theta*pi/180.0));
-      //  }
-    }
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: rnel/Kernel-Cylinder_f.cpp
===================================================================
--- Kernel/Kernel-Cylinder_f.cpp	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
+++ 	(revision )
@@ -1,113 +1,0 @@
-#ifdef __OPENCL_VERSION__
-# define USE_OPENCL
-#endif
-
-// If opencl is not available
-#ifndef USE_OPENCL
-#  include <math.h>
-#  define REAL(x) (x)
-#  define real double
-#  define global
-#  define local
-#  define kernel extern "C"
-#  include "NR_BessJ1.cpp"
-#endif
-
-#ifndef M_PI
-#  define M_PI REAL(3.141592653589793)
-#endif
-#define RADIANS REAL(0.017453292519943295)
-#define PI_E8 REAL(3.141592653589793e8)
-real f(real qx, real qy, real sub, real rr, real h, real scale, real cyl_theta, real cyl_phi);
-real f(real qx, real qy, real sub, real rr, real h, real scale, real cyl_theta, real cyl_phi)
-{
-    const real qq = sqrt(qx*qx+qy*qy);
-
-    const real theta = cyl_theta*RADIANS;
-    const real cos_val = cos(theta)*cos(cyl_phi*RADIANS)*(qx/qq) + sin(theta)*(qy/qq);
-
-    real alpha = acos(cos_val);
-    if(alpha == REAL(0.0)) {
-        alpha = REAL(1.0e-26);
-    }
-    const real sin_alpha = sin(alpha);
-    const real besarg = qq*rr*sin_alpha;
-    const real siarg = qq*h/2*cos(alpha);
-
-    real be, si;
-    if (besarg == REAL(0.0)) {
-        be = sin_alpha;
-    } else {
-        const real bj = NR_BessJ1(besarg)/besarg;
-        be = REAL(4.0)*bj*bj*sin_alpha;
-    }
-    if (siarg == REAL(0.0)) {
-        si = REAL(1.0);
-    } else {
-        si = sin(siarg)/siarg;
-    }
-    const real form = be*si*si/sin_alpha;
-    return PI_E8*scale*form*sub*sub*rr*rr*rr*rr*h*h;
-}
-
-kernel void CylinderKernel(global const real *qx, global const real *qy, global real *result,
-#ifdef USE_OPENCL
-    global real *loops_g,
-#else
-    const int Nq,
-#endif
-    local real *loops, const real cutoff,
-    const real scale, const real background,
-    const real sub,
-    const int Nradius, const int Nlength, const int Ntheta, const int Nphi)
-{
-#ifdef USE_OPENCL
-    // copy loops info to local memory
-    event_t e = async_work_group_copy(loops, loops_g, (Nradius+Nlength+Ntheta+Nphi)*2, 0);
-    wait_group_events(1, &e);
-
-    int i = get_global_id(0);
-    int count = get_global_size(0);
-
-    if(i < count)
-#else
-    #pragma omp parallel for
-    for (int i=0; i<Nq; i++)
-#endif
-    {
-        const real qxi=qx[i];
-        const real qyi=qy[i];
-        real ret=REAL(0.0), norm=REAL(0.0), norm_vol=REAL(0.0), vol=REAL(0.0);
-        for (int ri=0; ri < Nradius; ri++) {
-            const real rv = loops[2*ri];
-            const real rw = loops[2*ri+1];
-            for (int li=0; li < Nlength; li++) {
-                const real lv = loops[2*(li+Nradius)];
-                const real lw = loops[2*(li+Nradius)+1];
-                for (int thi=0; thi < Ntheta; thi++) {
-                    const real thv = loops[2*(thi+Nradius+Nlength)];
-                    const real thw = loops[2*(thi+Nradius+Nlength)+1];
-                    // #pragma unroll
-                    for (int phi=0; phi < Nphi; phi++) {
-                        const real phv = loops[2*(phi+Nradius+Nlength+Ntheta)];
-                        const real phw = loops[2*(phi+Nradius+Nlength+Ntheta)+1];
-
-                        const real weight = rw*lw*thw*phw;
-                        //ret += qxi + qyi + sub + rv + lv + weight + thv + phv;
-                        if (weight > cutoff) {
-                            ret += f(qxi, qyi, sub, rv, lv, weight, thv, phv);
-                            norm += weight;
-                            vol += rw*lw*rv*rv*lv;
-                            norm_vol += rw*lw;
-                        }
-                    }
-                }
-            }
-        }
-        //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: rnel/Kernel-Ellipse.cpp
===================================================================
--- Kernel/Kernel-Ellipse.cpp	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
+++ 	(revision )
@@ -1,29 +1,0 @@
-__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)
-{
-     int i = get_global_id(0);
-     if(i < length){
-         real ret = 0;
-         real q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
-         real pi = 4.0*atan(1.0);
-         real theta = axis_theta*pi/180.0;
-         real cos_val = cos(theta)*cos(axis_phi*pi/180.0)*(qx[i]/q) + sin(theta)*(qy[i]/q);
-
-         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*scale*(1.0e8);
-
-         _ptvalue[i] += radius_a_weight*radius_b_weight*axis_theta_weight*radius_a*axis_phi_weight*ret*pow(radius_b, 2);
-         //if(size > 1){
-          //  _ptvalue[i] *= fabs(cos(axis_theta*pi/180.0));
-         //}
-     }
-}
Index: rnel/Kernel-Ellipse_f.cpp
===================================================================
--- Kernel/Kernel-Ellipse_f.cpp	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
+++ 	(revision )
@@ -1,101 +1,0 @@
-#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];
-                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;
-                            vol += raw*rbw*rav*rbv*rbv;
-                            norm_vol += raw*rbw;
-                        }
-                    }
-                }
-            }
-        }
-        //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: rnel/Kernel-Lamellar.cpp
===================================================================
--- Kernel/Kernel-Lamellar.cpp	(revision 099e05357033bb53919d7a0213dc1bd7a0b54ebb)
+++ 	(revision )
@@ -1,11 +1,0 @@
-
-__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 int length)
-{
-    int i = get_global_id(0);
-    if(i < length)
-    {
-        real q = sqrt(qx[i]*qx[i]+qy[i]*qy[i]);
-        ret[i] = 2.0*4.0*atan(1.0)*scale*2.0*sub*(sub/q)/q*(1.0-cos(q*bi_thick))/(q*q)/bi_thick*1.0e8 ;
-    }
-}
Index: rnel/Kernel-OneDCylinder.cpp
===================================================================
--- Kernel/Kernel-OneDCylinder.cpp	(revision 099e05357033bb53919d7a0213dc1bd7a0b54ebb)
+++ 	(revision )
@@ -1,25 +1,0 @@
-__kernel void OneDCylKernel(__global const real *q, __global real *answer, const real sub, const real length,
-const real radius, const real scale, const real size, const real uplim, const real lolim)
-{
-    int k = get_global_id(0);
-    if(k < size)
-    {
-        real Pi = 4.0*atan(1.0);
-        real up = uplim*Pi/180;
-        real low = lolim*Pi/180;
-
-        real summ = 0.0; real zi = 0;
-        real halfheight = length/2.0;
-
-        for(int i=0; i<76; i++)
-        {
-            zi = (Gauss76Z(i)*(up-low) + up + low )/2.0;
-            summ += Gauss76Wt(i)*CylKernel(q[k], radius, halfheight, zi);
-        }
-
-        answer[k] = (up-low)/2.0*summ*sub*sub*1.0e8*scale*Pi*radius*radius*length;
-
-    }
-}
-
-
Index: rnel/Kernel-TriaxialEllipse.cpp
===================================================================
--- Kernel/Kernel-TriaxialEllipse.cpp	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ 	(revision )
@@ -1,43 +1,0 @@
-__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);
-        real answer=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 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: rnel/NR_BessJ1.cpp
===================================================================
--- Kernel/NR_BessJ1.cpp	(revision 723198c52509a7bfef1f123ffbabeebe978a7ad6)
+++ 	(revision )
@@ -1,23 +1,0 @@
-real NR_BessJ1(real x);
-real NR_BessJ1(real x)
-{
-    real xx=0; real y=0; real ans=0; real ans1=0; real ans2=0;
-    real ax = fabs(x);
-	if (ax < 8.0)
-	{
-		y=x*x;
-		ans1=x*(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 (x < 0.0) {ans*=-1;}
-	}
-	return(ans);
-}
Index: rnel/OneDCyl_Kfun.cpp
===================================================================
--- Kernel/OneDCyl_Kfun.cpp	(revision 099e05357033bb53919d7a0213dc1bd7a0b54ebb)
+++ 	(revision )
@@ -1,73 +1,0 @@
-real CylKernel(real qq, real rr, real h, real theta)
-{
-    // qq is the q-value for the calculation (1/A)
-    // rr is the radius of the cylinder (A)
-    // h is the HALF-LENGTH of the cylinder = L/2 (A)
-    real be=0; real si=0;
-    real besarg = qq*rr*sin(theta);
-    real siarg = qq*h*cos(theta);
-    real bj = NR_BessJ1(besarg);
-
-    real d1 = qq*rr*sin(theta);
-
-    if (besarg == 0.0){
-        be = sin(theta);
-    }else{
-        be = bj*bj*4.0*sin(theta) / (d1*d1);
-    }
-
-    if (siarg == 0.0){
-        si = 1.0;
-    }else{
-        si = sin(siarg)*sin(siarg)/(siarg*siarg);
-    }
-    return be*si;
-}
-real Gauss76Z(int i)
-{
-real array[76] = {.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) };
-
-return array[i];
-}
-real Gauss76Wt(int i)
-{
-real array[76] = {.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 };
-
-return array[i];
-}
Index: rnel/__init__.py
===================================================================
--- Kernel/__init__.py	(revision 473183c753f3b44a54aab3a16e00e358edf72be6)
+++ 	(revision )
@@ -1,1 +1,0 @@
-__author__ = 'HMP1'
Index: dels/__init__.py
===================================================================
--- Models/__init__.py	(revision 473183c753f3b44a54aab3a16e00e358edf72be6)
+++ 	(revision )
@@ -1,1 +1,0 @@
-__author__ = 'HMP1'
Index: dels/code_capcyl.py
===================================================================
--- Models/code_capcyl.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ 	(revision )
@@ -1,127 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-from math import sqrt, fabs, atan
-import pyopencl as cl
-
-from weights import GaussianDispersion
-from sasmodel import card, set_precision
-
-
-
-class GpuCapCylinder(object):
-    PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'len_cyl':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0,
-             'theta':0, 'phi':0}
-
-    PD_PARS = ['rad_cyl', 'len_cyl', 'rad_cap', 'theta', 'phi']
-
-    def __init__(self, qx, qy, dtype='float32'):
-
-        #create context, queue, and build program
-        ctx,_queue = card()
-        trala = open('Kernel/NR_BessJ1.cpp').read()+"\n"+open('Kernel/Capcyl_Kfun.cpp').read()+"\n"+open('Kernel/Kernel-CapCyl.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)
-        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()
-        self.res[:] = 0
-        cl.enqueue_copy(queue, self.res_b, self.res)
-
-        rad_cyl,len_cyl,rad_cap,theta,phi = \
-            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
-             for base in GpuCapCylinder.PD_PARS]
-
-        rad_cyl.value, rad_cyl.weight = rad_cyl.get_weights(pars['rad_cyl'], 0, 10000, True)
-        rad_cap.value, rad_cap.weight = rad_cap.get_weights(pars['rad_cap'], 0, 10000, True)
-        len_cyl.value, len_cyl.weight = len_cyl.get_weights(pars['len_cyl'], 0, 10000, True)
-        theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False)
-        phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False)
-
-        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
-        size = len(theta.weight)
-        sub = pars['sld_capcyl']-pars['sld_solv']
-        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-
-        for i in xrange(len(rad_cyl.weight)):
-            for m in xrange(len(rad_cap.weight)):
-                for j in xrange(len(len_cyl.weight)):
-
-                    hDist = -1.0*sqrt(fabs(rad_cap.value[m]*rad_cap.value[m]-rad_cyl.value[i]*rad_cyl.value[i]))
-                    vol_i = 4.0*atan(1.0)*rad_cyl.value[i]*rad_cyl.value[i]*len_cyl.value[j]+2.0*4.0*atan(1.0)/3.0\
-                                    *((rad_cap.value[m]-hDist)*(rad_cap.value[m]-hDist)*(2*rad_cap.value[m]+hDist))
-                    vol += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]*vol_i
-                    norm_vol += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]
-
-                    for k in xrange(len(theta.weight)):
-                        for l in xrange(len(phi.weight)):
-
-                            self.prg.CapCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                        real(vol_i), real(hDist), real(rad_cyl.value[i]), real(rad_cap.value[m]), real(len_cyl.value[j]),
-                                        real(theta.value[k]), real(phi.value[l]), real(sub), real(pars['scale']),
-                                        real(phi.weight[l]), real(theta.weight[k]), real(rad_cap.weight[m]),
-                                        real(rad_cyl.weight[i]), real(len_cyl.weight[j]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size))
-
-                            norm += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]*theta.weight[k]*phi.weight[l]
-
-
-
-        #if size > 1:
-         #   norm /= asin(1.0)
-        cl.enqueue_copy(queue, self.res, self.res_b)
-        sum += self.res
-        if vol != 0.0 and norm_vol != 0.0:
-            sum *= norm_vol/vol
-
-        return sum/norm + pars['background']
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: dels/code_coreshellcyl.py
===================================================================
--- Models/code_coreshellcyl.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ 	(revision )
@@ -1,77 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import pyopencl as cl
-
-from weights import GaussianDispersion
-from sasmodel import card, set_precision
-
-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('Kernel/NR_BessJ1.cpp').read()+"\n"+open('Kernel/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()
-        self.res[:] = 0
-        cl.enqueue_copy(queue, self.res_b, self.res)
-        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, 10000, True)
-        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
-        thickness.value, thickness.weight = thickness.get_weights(pars['thickness'], 0, 10000, 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)
-
-        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-        for r in xrange(len(radius.weight)):
-            for l in xrange(len(length.weight)):
-                for th in xrange(len(thickness.weight)):
-
-                    vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\
-                                   *(length.value[l]+2.0*thickness.value[th])
-                    norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th]
-
-                    for at in xrange(len(axis_theta.weight)):
-                        for p in xrange(len(axis_phi.weight)):
-
-                            self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                    real(axis_theta.value[at]), real(axis_phi.value[p]), real(thickness.value[th]),
-                                    real(length.value[l]), real(radius.value[r]), real(pars['scale']),
-                                    real(radius.weight[r]), real(length.weight[l]), real(thickness.weight[th]),
-                                    real(axis_theta.weight[at]), real(axis_phi.weight[p]), real(pars['core_sld']),
-                                    real(pars['shell_sld']), real(pars['solvent_sld']),np.uint32(size),
-                                    np.uint32(self.qx.size))
-
-                            norm += radius.weight[r]*length.weight[l]*thickness.weight[th]*axis_theta.weight[at]\
-                                    *axis_phi.weight[p]
-
-        #if size>1:
-         #   norm /= math.asin(1.0)
-        cl.enqueue_copy(queue, self.res, self.res_b)
-        sum = self.res
-        if vol != 0.0 and norm_vol != 0.0:
-            sum *= norm_vol/vol
-
-        return sum/norm + pars['background']
Index: dels/code_coreshellcyl_f.py
===================================================================
--- Models/code_coreshellcyl_f.py	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
+++ 	(revision )
@@ -1,92 +1,0 @@
-#!/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: dels/code_cylinder.py
===================================================================
--- Models/code_cylinder.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ 	(revision )
@@ -1,166 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import pyopencl as cl
-
-from weights import GaussianDispersion
-from sasmodel import card, set_precision, set_precision_1d, tic, toc
-
-
-class GpuCylinder(object):
-    PARS = {
-        'scale':1,'radius':1,'length':1,'sldCyl':1e-6,'sldSolv':0,'background':0,
-        'cyl_theta':0,'cyl_phi':0,
-    }
-    PD_PARS = ['radius', 'length', 'cyl_theta', 'cyl_phi']
-
-    def __init__(self, qx, qy, dtype='float32'):
-
-        #create context, queue, and build program
-        ctx,_queue = card()
-        src, qx, qy = set_precision(open('Kernel/NR_BessJ1.cpp').read()+"\n"+open('Kernel/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, cl.mem_flags.READ_WRITE, self.qx.nbytes)
-        self.res = np.empty_like(self.qx)
-
-    def eval(self, pars):
-
-        tic()
-        _ctx,queue = card()
-        self.res[:] = 0
-        cl.enqueue_copy(queue, self.res_b, self.res)
-        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, 10000, True)
-        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
-        cyl_theta.value, cyl_theta.weight = cyl_theta.get_weights(pars['cyl_theta'], -np.inf, np.inf, False)
-        cyl_phi.value, cyl_phi.weight = cyl_phi.get_weights(pars['cyl_phi'], -np.inf, np.inf, 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)):
-
-                vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j]
-                norm_vol += radius.weight[i]*length.weight[j]
-
-                for k in xrange(len(cyl_theta.weight)):
-                    for l in xrange(len(cyl_phi.weight)):
-                        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))
-
-                        norm += radius.weight[i]*length.weight[j]*cyl_theta.weight[k]*cyl_phi.weight[l]
-
-
-       # if size > 1:
-        #    norm /= math.asin(1.0)
-        cl.enqueue_copy(queue, self.res, self.res_b)
-        sum = self.res
-        if vol != 0.0 and norm_vol != 0.0:
-            sum *= norm_vol/vol
-
-        print toc()*1000, self.qx.shape[0]
-        return sum/norm+pars['background']
-
-class OneDGpuCylinder(object):
-    PARS = {
-        'scale':1,'radius':1,'length':1,'sldCyl':1e-6,'sldSolv':0,'background':0,
-        'bolim':0, 'uplim':90
-    }
-    PD_PARS = ['radius', 'length']
-
-    def __init__(self, q, dtype='float32'):
-
-        #create context, queue, and build program
-        ctx,_queue = card()
-        trala = open('Kernel/NR_BessJ1.cpp').read()+"\n"+open('Kernel/OneDCyl_Kfun.cpp').read()+"\n"+open('Kernel/Kernel-OneDCylinder.cpp').read()
-        src, self.q = set_precision_1d(trala, q, dtype=dtype)
-        self.prg = cl.Program(ctx, src).build()
-
-        #buffers
-        mf = cl.mem_flags
-        self.q_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.q)
-        self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, q.nbytes)
-        self.res = np.empty_like(self.q)
-
-    def eval(self, pars):
-
-        _ctx,queue = card()
-        radius, length = \
-            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
-             for base in OneDGpuCylinder.PD_PARS]
-
-        #Get the weights for each
-        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True)
-        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
-
-        #Perform the computation, with all weight points
-        sum, norm, vol = 0.0, 0.0, 0.0,
-        sub = pars['sldCyl'] - pars['sldSolv']
-
-        real = np.float32 if self.q.dtype == np.dtype('float32') else np.float64
-        #Loop over radius, length, theta, phi weight points
-        for r in xrange(len(radius.weight)):
-            for l in xrange(len(length.weight)):
-                        self.prg.OneDCylKernel(queue, self.q.shape, None, self.q_b, self.res_b, real(sub),
-                                           real(length.value[l]), real(radius.value[r]), real(pars['scale']),
-                                           np.uint32(self.q.size), real(pars['uplim']), real(pars['bolim']))
-                        cl.enqueue_copy(queue, self.res, self.res_b)
-                        sum += radius.weight[r]*length.weight[l]*self.res*pow(radius.value[r],2)*length.value[l]
-                        vol += radius.weight[r]*length.weight[l] *pow(radius.value[r],2)*length.value[l]
-                        norm += radius.weight[r]*length.weight[l]
-
-        if vol != 0.0 and norm != 0.0:
-            sum *= norm/vol
-
-        return sum/norm + pars['background']
-
-       
-
-        
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: dels/code_cylinder_f.py
===================================================================
--- Models/code_cylinder_f.py	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
+++ 	(revision )
@@ -1,149 +1,0 @@
-#!/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 GpuCylinder(object):
-    PARS = {
-        'scale':1,'radius':1,'length':1,'sldCyl':1e-6,'sldSolv':0,'background':0,
-        'cyl_theta':0,'cyl_phi':0,
-    }
-    PD_PARS = ['radius', 'length', 'cyl_theta', 'cyl_phi']
-
-    def __init__(self, qx, qy, dtype='float32', cutoff=1e-5):
-
-        #create context, queue, and build program
-        ctx,_queue = card()
-        bessel = open('Kernel/NR_BessJ1.cpp').read()
-        kernel = open('Kernel/Kernel-Cylinder_f.cpp').read()
-        src, qx, qy = set_precision("\n".join((bessel,kernel)), qx, qy, dtype=dtype)
-        self.prg = cl.Program(ctx, src).build()
-        self.qx, self.qy = qx, qy
-        self.cutoff = cutoff
-
-        #buffers
-        self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
-        self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
-        self.res_b = cl.Buffer(ctx, cl.mem_flags.READ_WRITE, self.qx.nbytes)
-        self.res = np.empty_like(self.qx)
-
-    def eval(self, pars):
-        tic()
-
-        ctx,queue = card()
-        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-        loops, loop_lengths = make_loops(pars, dtype=self.qx.dtype)
-        loops_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops)
-        loops_l = cl.LocalMemory(len(loops.data))
-
-        self.prg.CylinderKernel(queue, self.qx.shape, None,
-            self.qx_b, self.qy_b, self.res_b,
-            loops_b, loops_l, real(self.cutoff),
-            real(pars['scale']), real(pars['background']),
-            real(pars['sldCyl']-pars['sldSolv']),
-            *[np.uint32(pn) for pn in loop_lengths])
-
-        cl.enqueue_copy(queue, self.res, self.res_b)
-        print toc()*1000, self.qx.shape[0]
-
-        return self.res
-
-class CpuCylinder(GpuCylinder): #inherit parameters only
-    def __init__(self, qx, qy, dtype='float32', cutoff=1e-5):
-        self.qx, self.qy = [np.ascontiguousarray(v,'d') for v in qx,qy]
-        self.cutoff = cutoff
-        self.res = np.empty_like(self.qx)
-        self.dll = ctypes.CDLL('Kernel/cylinder.so')
-        self.fn = self.dll['CylinderKernel']
-        self.fn.argtypes = [
-            c_void_p, c_void_p, c_void_p, c_int,
-            c_void_p, c_double, c_double, c_double, c_double,
-            c_int, c_int, c_int, c_int
-        ]
-    def eval(self, pars):
-        loops, loop_lengths = make_loops(pars, dtype=self.qx.dtype)
-        self.fn(self.qx.ctypes.data, self.qy.ctypes.data, self.res.ctypes.data, len(self.qx),
-            loops.ctypes.data, self.cutoff, pars['scale'], pars['background'],
-            pars['sldCyl']-pars['sldSolv'], *loop_lengths)
-
-        return self.res
-
-def make_loops(pars, dtype='double'):
-    # 0.2 ms on sparkle to form the final loops
-    radius, length, theta, phi = \
-       [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
-         for base in GpuCylinder.PD_PARS]
-    parts = [
-        radius.get_weights(pars['radius'], 0, 10000, True),
-        length.get_weights(pars['length'], 0, 10000, True),
-        theta.get_weights(pars['cyl_theta'], -np.inf, np.inf, False),
-        phi.get_weights(pars['cyl_phi'], -np.inf, np.inf, False),
-        ]
-    # Make sure that weights are normalized to peaks at 1 so that
-    # the tolerance term can be used properly on truncated distributions
-    loops = np.hstack((v,w/w.max()) for v,w in parts)
-    #loops = np.hstack(parts)
-    loops = np.ascontiguousarray(loops.T, dtype).flatten()
-    return loops, [len(p[0]) for p in parts]
-
-class OneDGpuCylinder(object):
-    PARS = {
-        'scale':1,'radius':1,'length':1,'sldCyl':1e-6,'sldSolv':0,'background':0,
-        'bolim':0, 'uplim':90
-    }
-    PD_PARS = ['radius', 'length']
-
-    def __init__(self, q, dtype='float32'):
-
-        #create context, queue, and build program
-        ctx,_queue = card()
-        trala = open('Kernel/NR_BessJ1.cpp').read()+"\n"+open('Kernel/OneDCyl_Kfun.cpp').read()+"\n"+open('Kernel/Kernel-OneDCylinder.cpp').read()
-        src, self.q = set_precision_1d(trala, q, dtype=dtype)
-        self.prg = cl.Program(ctx, src).build()
-
-        #buffers
-        mf = cl.mem_flags
-        self.q_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.q)
-        self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, q.nbytes)
-        self.res = np.empty_like(self.q)
-
-    def eval(self, pars):
-
-        _ctx,queue = card()
-        radius, length = \
-            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
-             for base in OneDGpuCylinder.PD_PARS]
-
-        #Get the weights for each
-        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True)
-        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
-
-        #Perform the computation, with all weight points
-        sum, norm, vol = 0.0, 0.0, 0.0,
-        sub = pars['sldCyl'] - pars['sldSolv']
-
-        real = np.float32 if self.q.dtype == np.dtype('float32') else np.float64
-        #Loop over radius, length, theta, phi weight points
-        for r in xrange(len(radius.weight)):
-            for l in xrange(len(length.weight)):
-                        self.prg.OneDCylKernel(queue, self.q.shape, None, self.q_b, self.res_b, real(sub),
-                                           real(length.value[l]), real(radius.value[r]), real(pars['scale']),
-                                           np.uint32(self.q.size), real(pars['uplim']), real(pars['bolim']))
-                        cl.enqueue_copy(queue, self.res, self.res_b)
-                        sum += radius.weight[r]*length.weight[l]*self.res*pow(radius.value[r],2)*length.value[l]
-                        vol += radius.weight[r]*length.weight[l] *pow(radius.value[r],2)*length.value[l]
-                        norm += radius.weight[r]*length.weight[l]
-
-        if vol != 0.0 and norm != 0.0:
-            sum *= norm/vol
-
-        return sum/norm + pars['background']
Index: dels/code_ellipse.py
===================================================================
--- Models/code_ellipse.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ 	(revision )
@@ -1,119 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import pyopencl as cl
-
-from weights import GaussianDispersion
-from sasmodel import card, set_precision
-
-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'):
-
-        ctx,_queue = card()
-        src, qx, qy = set_precision(open('Kernel/Kernel-Ellipse.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):
-    #b_n = radius_b # want, a_n = radius_a # want, etc
-        _ctx,queue = card()
-        self.res[:] = 0
-        cl.enqueue_copy(queue, self.res_b, self.res)
-        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, 10000, True)
-        radius_b.value, radius_b.weight = radius_b.get_weights(pars['radius_b'], 0, 10000, True)
-        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False)
-        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False)
-
-        #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']
-        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-
-        #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
-                vol += radius_a.weight[i]*radius_b.weight[j]*pow(radius_b.value[j], 2)*radius_a.value[i]
-                norm_vol += radius_a.weight[i]*radius_b.weight[j]
-
-                for k in xrange(len(axis_theta.weight)):
-                    #Average over phi distribution
-                    for l in xrange(len(axis_phi.weight)):
-                        #call the kernel
-                        self.prg.EllipsoidKernel(queue, self.qx.shape, None, real(radius_a.weight[i]),
-                                        real(radius_b.weight[j]), real(axis_theta.weight[k]),
-                                        real(axis_phi.weight[l]), real(pars['scale']), real(radius_a.value[i]),
-                                        real(radius_b.value[j]), real(sub), real(axis_theta.value[k]),
-                                        real(axis_phi.value[l]), self.qx_b, self.qy_b, self.res_b,
-                                        np.uint32(self.qx.size), np.uint32(len(axis_theta.weight)))
-
-                        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)
-        cl.enqueue_copy(queue, self.res, self.res_b)
-        sum += self.res
-        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()
-
-    #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: dels/code_ellipse_f.py
===================================================================
--- Models/code_ellipse_f.py	(revision ae7d639faa3ed3b76696623477e23d48d49ddeff)
+++ 	(revision )
@@ -1,89 +1,0 @@
-#!/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: dels/code_lamellar.py
===================================================================
--- Models/code_lamellar.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ 	(revision )
@@ -1,77 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import pyopencl as cl
-from Models.weights import GaussianDispersion
-from sasmodel import set_precision
-
-
-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/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, 10000, 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']
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: dels/code_triaxialellipse.py
===================================================================
--- Models/code_triaxialellipse.py	(revision 1726b21cbb8c6f89faab2157203ecaa4b5980db8)
+++ 	(revision )
@@ -1,118 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import pyopencl as cl
-
-from weights import GaussianDispersion
-from sasmodel import card, set_precision
-
-class GpuTriEllipse:
-    PARS = {'scale':1, 'semi_axisA':35, 'semi_axisB':100, 'semi_axisC':400, 'sldEll':1e-6, 'sldSolv':6.3e-6, 'background':0,
-            'axis_theta':0, 'axis_phi':0, 'axis_psi':0}
-
-    PD_PARS = ['semi_axisA', 'semi_axisB', 'semi_axisC', 'axis_theta', 'axis_phi', 'axis_psi']
-
-    def __init__(self, qx, qy, dtype='float32'):
-        ctx,_queue = card()
-        src, qx, qy = set_precision(open('Kernel/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()
-        self.res[:] = 0
-        cl.enqueue_copy(queue, self.res_b, self.res)
-        semi_axisA, semi_axisB, semi_axisC, axis_theta, axis_phi, axis_psi = \
-            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
-             for base in GpuTriEllipse.PD_PARS]
-
-        semi_axisA.value, semi_axisA.weight = semi_axisA.get_weights(pars['semi_axisA'], 0, 10000, True)
-        semi_axisB.value, semi_axisB.weight = semi_axisB.get_weights(pars['semi_axisB'], 0, 10000, True)
-        semi_axisC.value, semi_axisC.weight = semi_axisC.get_weights(pars['semi_axisC'], 0, 10000, True)
-        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False)
-        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False)
-        axis_psi.value, axis_psi.weight = axis_psi.get_weights(pars['axis_psi'], -90, 180, False)
-
-        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
-        size = len(axis_theta.weight)
-        sub = pars['sldEll'] - pars['sldSolv']
-
-        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-        for a in xrange(len(semi_axisA.weight)):
-            for b in xrange(len(semi_axisB.weight)):
-                for c in xrange(len(semi_axisC.weight)):
-
-                    vol += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]*semi_axisA.value[a]*semi_axisB.value[b]*semi_axisC.value[c]
-                    norm_vol += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]
-
-                    for t in xrange(len(axis_theta.weight)):
-                        for i in xrange(len(axis_phi.weight)):
-                            for s in xrange(len(axis_psi.weight)):
-                                self.prg.TriaxialEllipseKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                            real(sub), real(pars['scale']), real(semi_axisA.value[a]), real(semi_axisB.value[b]),
-                                            real(semi_axisC.value[c]), real(axis_phi.value[i]), real(axis_theta.value[t]), real(axis_psi.value[s]),
-                                            real(semi_axisA.weight[a]), real(semi_axisB.weight[b]), real(semi_axisC.weight[c]), real(axis_psi.weight[s]),
-                                            real(axis_phi.weight[i]), real(axis_theta.weight[t]), np.uint32(self.qx.size), np.uint32(size))
-
-                                norm += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]*axis_theta.weight[t]*axis_phi.weight[i]*axis_psi.weight[s]
-
-
-      #  if size > 1:
-       #     norm /= asin(1.0)
-        cl.enqueue_copy(queue, self.res, self.res_b)
-        sum = self.res
-        if vol != 0.0 and norm_vol != 0.0:
-            sum *= norm_vol/vol
-
-        return sum/norm + pars['background']
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: dels/weights.py
===================================================================
--- Models/weights.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ 	(revision )
@@ -1,26 +1,0 @@
-import numpy as np
-
-class GaussianDispersion(object):
-    def __init__(self, npts=35, width=0, nsigmas=3): #number want, percent deviation, #standard deviations from mean
-        self.type = 'gaussian'
-        self.npts = npts
-        self.width = width
-        self.nsigmas = nsigmas
-
-    def get_pars(self):
-        return self.__dict__
-
-    def get_weights(self, center, min, max, relative):
-        """ *center* is the center of the distribution
-        *min*,*max* are the min, max allowed values
-        *relative* is True if the width is relative to the center instead of absolute
-        For polydispersity use relative.  For orientation parameters use absolute."""
-        npts, width, nsigmas = self.npts, self.width, self.nsigmas
-        sigma = width * center if relative else width
-        if sigma == 0:
-            return np.array([center],'d'), np.array([1.], 'd')
-        x = center + np.linspace(-nsigmas * sigma, +nsigmas * sigma, npts)
-        x = x[(x >= min) & (x <= max)]
-        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
-        return x, px
-
Index: sters/TEST-Kernel-Ellipse.cpp
===================================================================
--- Testers/TEST-Kernel-Ellipse.cpp	(revision 099e05357033bb53919d7a0213dc1bd7a0b54ebb)
+++ 	(revision )
@@ -1,81 +1,0 @@
-__kernel void EllipsoidKernel(__global const real *qx, __global const real *qy, __global const int *place, __global const real *array,
-__global real *final, const real scale, const real sub, const int length, const int size)//, __local real *poly)
-{
-    //__local real rada, radaw, radb, radbw, th, thw, ph, phw;
-
-    __local real rada;
-    int uno=0; int dos=1; int count = 0;
-    int l = get_local_id(0);
-    for(int j = 0; j < 4; j++)
-    {
-        for(int i = place[uno]; i < place[dos]; i++)
-        {
-            rada[count] = array[i]
-            radaw[count] = array[]
-        }
-    uno+=2;
-    dos+=2;
-    count = 0;
-    }
-
-
-    int i = get_global_id(0);
-    if(i < length){
-    //do local mem?
-    real sum=0.0; real vol=0.0; real norm_vol=0.0; real norm=0.0;
-
-    for(int a = 0; a < place[0]; a++) //radius_a values
-    {
-        real radius_a = array[a];
-        real radius_a_weight = array[a+place[0]];
-        for(int b = place[1]; b < place[2]; b++) //radius_b values
-        {
-            real radius_b = array[b];
-            real radius_b_weight = array[b+place[2]];
-            for(int t = place[3]; t < place[4]; t++) //Axis_theta values
-            {
-                real axis_theta = array[t];
-                real axis_theta_weight = array[t+place[4]];
-                for(int p = place[5]; p < place[6]; p++) //axis_phi values
-                {
-                    real axis_phi = array[p];
-                    real axis_phi_weight = array[p+place[6]];
-
-                    real ret = 0;
-                    real q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
-                    real pi = 4.0*atan(1.0);
-                    real theta = axis_theta*pi/180.0;
-                    real cos_val = cos(theta)*cos(axis_phi*pi/180.0)*(qx[i]/q) + sin(theta)*(qy[i]/q);
-
-                    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*scale*(1.0e8);
-                    real _ptvalue = radius_a_weight*radius_b_weight*axis_theta_weight*radius_a*axis_phi_weight*ret*pow(radius_b, 2);
-                    if(size > 1){
-                        _ptvalue *= fabs(cos(axis_theta*pi/180.0));
-                    }
-
-                    sum += _ptvalue;
-                    vol += radius_a_weight*radius_b_weight*pow(radius_b, 2)*radius_a;
-                    norm_vol += radius_a_weight*radius_b_weight;
-                    norm += radius_a_weight*radius_b_weight*axis_theta_weight*axis_phi_weight;
-                }
-            }
-        }
-    }
-
-    if(size > 1){
-        norm /= asin(1.0);
-    }
-    if(vol != 0.0 && norm_vol != 0.0){
-        sum *= norm_vol/vol;
-    }
-
-    final[i] = sum/norm;
-}
-}
Index: sters/__init__.py
===================================================================
--- Testers/__init__.py	(revision 473183c753f3b44a54aab3a16e00e358edf72be6)
+++ 	(revision )
@@ -1,1 +1,0 @@
-__author__ = 'HMP1'
Index: sters/load.py
===================================================================
--- Testers/load.py	(revision 473183c753f3b44a54aab3a16e00e358edf72be6)
+++ 	(revision )
@@ -1,70 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-from sans.dataloader.loader import Loader
-from sans.dataloader.manipulations import Ringcut
-from Kernel import fit
-
-
-def cylinder(data): ##################################################################################
-    from sans.models.CylinderModel import CylinderModel
-    model1 = CylinderModel()
-    model1.setParam("scale", 10.0)
-    model1.setParam("radius",18)
-    model1.setParam("length", 397)
-    model1.setParam("sldCyl",3e-006 )
-    model1.setParam("sldSolv",0.0 )
-    model1.setParam("background", 0.0)
-
-    # Dispersion parameters
-    model1.dispersion['radius']['width'] = 0.25
-    model1.dispersion['radius']['npts'] = 50
-
-    theory = model1.evalDistribution([data.qx_data, data.qy_data])
-    return theory
-
-def load_data(filename):
-    loader = Loader()
-    f = loader.load(filename)
-    return f
-
-def set_beam_stop(data, radius):
-    data.mask = Ringcut(0, radius)(data)
-
-def plot_data(data):
-    from numpy.ma import masked_array
-    import matplotlib.pyplot as plt
-    img = masked_array(data.data, data.mask)
-    xmin, xmax = min(data.qx_data), max(data.qx_data)
-    ymin, ymax = min(data.qy_data), max(data.qy_data)
-    plt.imshow(img.reshape(128,128),
-               interpolation='nearest', aspect=1, origin='upper',
-               extent=[xmin, xmax, ymin, ymax])
-
-def demo():
-
-    data = load_data('NOV07090.DAT')
-    """
-    print data
-    print type(data)
-    set_beam_stop(data, 0.004)
-    plot_data(data)
-    import matplotlib.pyplot as plt; plt.show()
-    """
-    import matplotlib.pyplot as plt
-    import numpy as np
-    sasviewcyl = cylinder(data)
-    gpucyl = fit.Fit(data)
-    diff = gpucyl - sasviewcyl
-    np.linalg.norm(diff, order=None)
-    np.linalg.norm(diff, order=np.inf)
-
-    plot_data(data, diff)
-    plt.show()
-
-
-
-
-if __name__ == "__main__":
-    demo()
-
Index: sters/test.py
===================================================================
--- Testers/test.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ 	(revision )
@@ -1,19 +1,0 @@
-impor
-
-
-
-
-import numpy as np
-
-
-x = np.array([0, 1, 2, 3])
-y = 2*x + 1
-
-print y
-sum = 0
-
-for i in xrange(x.size-1):
-    sum += (x[i+1]-x[i])*(y[i]+y[i+1])/2
-
-print sum
-
Index: sters/test_code_ellipse.py
===================================================================
--- Testers/test_code_ellipse.py	(revision 099e05357033bb53919d7a0213dc1bd7a0b54ebb)
+++ 	(revision )
@@ -1,89 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import numpy as np
-import pyopencl as cl
-
-from weights import GaussianDispersion
-from Models.sasmodel import card
-import hi
-
-
-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 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'):
-
-        ctx,_queue = card()
-        src, qx, qy = set_precision(open('TEST-Kernel-Ellipse.cpp').read(), qx, qy, dtype=dtype)
-        self.prg = cl.Program(ctx, src).build()
-        self.qx, self.qy = qx, qy
-        place = np.ascontiguousarray(hi.place, dtype=int)
-        #buffers
-        mf = cl.mem_flags
-        self.place_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=place)
-        self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
-        self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
-        self.res_b = cl.Buffer(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
-        ctx,queue = card()
-        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']
-        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-
-        x = [radius_a.value, radius_a.weight, radius_b.value, radius_b.weight, axis_theta.value,
-                 axis_theta.weight, axis_phi.value, axis_phi.weight]
-        array = np.hstack(x)
-
-        array_b = cl.Buffer(ctx, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=array)
-
-        self.prg.EllipsoidKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.place_b, array_b, self.res_b,
-                                 real(pars['scale']), real(sub), np.uint32(self.qx.size), np.uint32(len(axis_theta.weight)))
-        #copy result back from buffer
-        cl.enqueue_copy(queue, self.res, self.res_b)
-
-
-        a = open("answer.txt", "w")
-        for x in xrange(len(self.res)):
-            a.write(str(self.res))
-            a.write("\n")
-
-
-        return self.res+pars['background']
-
-
-
-
-
Index: mpare-new.py
===================================================================
--- compare-new.py	(revision f4cf580e329add99780558498dbc8a8ab076ad24)
+++ 	(revision )
@@ -1,405 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import sys
-import math
-
-import numpy as np
-
-from sasmodels.bumps_model import BumpsModel, plot_data, tic
-from sasmodels import gpu, dll
-from sasmodels.convert import revert_model
-
-
-def sasview_model(modelname, **pars):
-    """
-    Load a sasview model given the model name.
-    """
-    # convert model parameters from sasmodel form to sasview form
-    #print "old",sorted(pars.items())
-    modelname, pars = revert_model(modelname, pars)
-    #print "new",sorted(pars.items())
-    sans = __import__('sans.models.'+modelname)
-    ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None)
-    if ModelClass is None:
-        raise ValueError("could not find model %r in sans.models"%modelname)
-    model = ModelClass()
-
-    for k,v in pars.items():
-        if k.endswith("_pd"):
-            model.dispersion[k[:-3]]['width'] = v
-        elif k.endswith("_pd_n"):
-            model.dispersion[k[:-5]]['npts'] = v
-        elif k.endswith("_pd_nsigma"):
-            model.dispersion[k[:-10]]['nsigmas'] = v
-        elif k.endswith("_pd_type"):
-            model.dispersion[k[:-8]]['type'] = v
-        else:
-            model.setParam(k, v)
-    return model
-
-def load_opencl(modelname, dtype='single'):
-    sasmodels = __import__('sasmodels.models.'+modelname)
-    module = getattr(sasmodels.models, modelname, None)
-    kernel = gpu.load_model(module, dtype=dtype)
-    return kernel
-
-def load_ctypes(modelname, dtype='single'):
-    sasmodels = __import__('sasmodels.models.'+modelname)
-    module = getattr(sasmodels.models, modelname, None)
-    kernel = dll.load_model(module, dtype=dtype)
-    return kernel
-
-def randomize(p, v):
-    """
-    Randomizing parameter.
-
-    Guess the parameter type from name.
-    """
-    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')):
-        return v
-    elif any(s in p for s in ('theta','phi','psi')):
-        # orientation in [-180,180], orientation pd in [0,45]
-        if p.endswith('_pd'):
-            return 45*np.random.rand()
-        else:
-            return 360*np.random.rand() - 180
-    elif 'sld' in p:
-        # sld in in [-0.5,10]
-        return 10.5*np.random.rand() - 0.5
-    elif p.endswith('_pd'):
-        # length pd in [0,1]
-        return np.random.rand()
-    else:
-        # length, scale, background in [0,200]
-        return 200*np.random.rand()
-
-def parlist(pars):
-    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
-
-def suppress_pd(pars):
-    """
-    Suppress theta_pd for now until the normalization is resolved.
-
-    May also suppress complete polydispersity of the model to test
-    models more quickly.
-    """
-    for p in pars:
-        if p.endswith("_pd"): pars[p] = 0
-
-def compare(name, pars, Ncpu, Ngpu, opts, set_pars):
-
-    # Sort out data
-    qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05)
-    if "-1d" in opts:
-        from sasmodels.bumps_model import empty_data1D
-        qmax = math.log10(qmax)
-        data = empty_data1D(np.logspace(qmax-3, qmax, 128))
-        index = slice(None, None)
-    else:
-        from sasmodels.bumps_model import empty_data2D, set_beam_stop
-        data = empty_data2D(np.linspace(-qmax, qmax, 128))
-        set_beam_stop(data, 0.004)
-        index = ~data.mask
-    is2D = hasattr(data, 'qx_data')
-
-    # modelling accuracy is determined by dtype and cutoff
-    dtype = 'double' if '-double' in opts else 'single'
-    cutoff_opts = [s for s in opts if s.startswith('-cutoff')]
-    cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5
-
-    # randomize parameters
-    random_opts = [s for s in opts if s.startswith('-random')]
-    if random_opts:
-        if '=' in random_opts[0]:
-            seed = int(random_opts[0].split('=')[1])
-        else:
-            seed = int(np.random.rand()*10000)
-        np.random.seed(seed)
-        print "Randomize using -random=%i"%seed
-        # Note: the sort guarantees order of calls to random number generator
-        pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items()))
-        # The capped cylinder model has a constraint on its parameters
-        if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
-            pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
-    pars.update(set_pars)
-
-    # parameter selection
-    if '-mono' in opts:
-        suppress_pd(pars)
-    if '-pars' in opts:
-        print "pars",parlist(pars)
-
-    # OpenCl calculation
-    if Ngpu > 0:
-        try:
-            model = load_opencl(name, dtype=dtype)
-        except Exception,exc:
-            print exc
-            print "... trying again with single precision"
-            model = load_opencl(name, dtype='single')
-        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
-        toc = tic()
-        for _ in range(Ngpu):
-            #pars['scale'] = np.random.rand()
-            problem.update()
-            gpu = problem.theory()
-        gpu_time = toc()*1000./Ngpu
-        print "opencl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[index]))
-        #print max(gpu), min(gpu)
-
-    # ctypes/sasview calculation
-    if Ncpu > 0 and "-ctypes" in opts:
-        model = load_ctypes(name, dtype=dtype)
-        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
-        toc = tic()
-        for _ in range(Ncpu):
-            problem.update()
-            cpu = problem.theory()
-        cpu_time = toc()*1000./Ncpu
-        comp = "ctypes"
-        print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
-    elif Ncpu > 0:
-        model = sasview_model(name, **pars)
-        toc = tic()
-        for _ in range(Ncpu):
-            if is2D:
-                cpu = model.evalDistribution([data.qx_data, data.qy_data])
-            else:
-                cpu = model.evalDistribution(data.x)
-        cpu_time = toc()*1000./Ncpu
-        comp = "sasview"
-        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
-
-    # Compare, but only if computing both forms
-    if Ngpu > 0 and Ncpu > 0:
-        #print "speedup %.2g"%(cpu_time/gpu_time)
-        #print "max |gpu/cpu|", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu))
-        #cpu *= max(gpu/cpu)
-        resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu)
-        resid[index] = (gpu - cpu)[index]
-        relerr[index] = resid[index]/cpu[index]
-        print "max(|ocl-%s|)"%comp, max(abs(resid[index]))
-        print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[index]))
-
-    # Plot if requested
-    if '-noplot' in opts: return
-    import matplotlib.pyplot as plt
-    if Ncpu > 0:
-        if Ngpu > 0: plt.subplot(131)
-        plot_data(data, cpu, scale='log')
-        plt.title("%s t=%.1f ms"%(comp,cpu_time))
-    if Ngpu > 0:
-        if Ncpu > 0: plt.subplot(132)
-        plot_data(data, gpu, scale='log')
-        plt.title("opencl t=%.1f ms"%gpu_time)
-    if Ncpu > 0 and Ngpu > 0:
-        plt.subplot(133)
-        err = resid if '-abs' in opts else relerr
-        errstr = "abs err" if '-abs' in opts else "rel err"
-        #err,errstr = gpu/cpu,"ratio"
-        plot_data(data, err, scale='linear')
-        plt.title("max %s = %.3g"%(errstr, max(abs(err[index]))))
-        if is2D: plt.colorbar()
-    plt.show()
-
-# ===========================================================================
-#
-USAGE="""
-usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val]
-
-Compare the speed and value for a model between the SasView original and the
-OpenCL rewrite.
-
-model is the name of the model to compare (see below).
-Nopencl is the number of times to run the OpenCL model (default=5)
-Nsasview is the number of times to run the Sasview model (default=1)
-
-Options (* for default):
-
-    -plot*/-noplot plots or suppress the plot of the model
-    -single/-double uses double precision for comparison
-    -lowq/-midq/-highq use q values up to 0.05, 0.2 or 1.0
-    -2d/-1d uses 1d or 2d random data
-    -preset/-random[=seed] preset or random parameters
-    -poly/-mono force monodisperse/polydisperse
-    -sasview/-ctypes whether cpu is tested using sasview or ctypes
-    -cutoff=1e-5/value cutoff for including a point in polydispersity
-    -nopars/-pars prints the parameter set or not
-    -rel/-abs plot relative or absolute error
-
-Key=value pairs allow you to set specific values to any of the model
-parameters.
-
-Available models:
-
-    %s
-"""
-
-VALID_OPTIONS = [
-    'plot','noplot',
-    'single','double',
-    'lowq','midq','highq',
-    '2d','1d',
-    'preset','random',
-    'poly','mono',
-    'sasview','ctypes',
-    'nopars','pars',
-    'rel','abs',
-    ]
-
-def main():
-    opts = [arg for arg in sys.argv[1:] if arg.startswith('-')]
-    args = [arg for arg in sys.argv[1:] if not arg.startswith('-')]
-    models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' '))
-                           for k,v in sorted(MODELS.items()))
-    if len(args) == 0:
-        print(USAGE%models)
-        sys.exit(1)
-    if args[0] not in MODELS:
-        print "Model %r not available. Use one of:\n    %s"%(args[0],models)
-        sys.exit(1)
-
-    valid_opts = set(VALID_OPTIONS)
-    invalid = [o[1:] for o in opts
-               if o[1:] not in valid_opts
-                  and not o.startswith('-cutoff=')
-                  and not o.startswith('-random=')]
-    if invalid:
-        print "Invalid options: %s"%(", ".join(invalid))
-        sys.exit(1)
-
-    name, pars = MODELS[args[0]]()
-    Nopencl = int(args[1]) if len(args) > 1 else 5
-    Nsasview = int(args[2]) if len(args) > 3 else 1
-
-    # Fill in default polydispersity parameters
-    pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd'))
-    for p in pds:
-        if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3
-        if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian"
-
-    # Fill in parameters given on the command line
-    set_pars = {}
-    for arg in args[3:]:
-        k,v = arg.split('=')
-        if k not in pars:
-            # extract base name without distribution
-            s = set(p.split('_pd')[0] for p in pars)
-            print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s)))
-            sys.exit(1)
-        set_pars[k] = float(v) if not v.endswith('type') else v
-
-    compare(name, pars, Nsasview, Nopencl, opts, set_pars)
-
-# ===========================================================================
-#
-
-MODELS = {}
-def model(name):
-    def gather_function(fn):
-        MODELS[name] = fn
-        return fn
-    return gather_function
-
-
-@model('cyl')
-def cylinder():
-    pars = dict(
-        scale=1, background=0,
-        sld=6, solvent_sld=1,
-        #radius=5, length=20,
-        radius=260, length=290,
-        theta=30, phi=0,
-        radius_pd=.2, radius_pd_n=1,
-        length_pd=.2,length_pd_n=1,
-        theta_pd=15, theta_pd_n=45,
-        phi_pd=15, phi_pd_n=1,
-        )
-    return 'cylinder', pars
-
-@model('capcyl')
-def capped_cylinder():
-    pars = dict(
-        scale=1, background=0,
-        sld=6, solvent_sld=1,
-        radius=260, cap_radius=290, length=290,
-        theta=30, phi=15,
-        radius_pd=.2, radius_pd_n=1,
-        cap_radius_pd=.2, cap_radius_pd_n=1,
-        length_pd=.2, length_pd_n=1,
-        theta_pd=15, theta_pd_n=45,
-        phi_pd=15, phi_pd_n=1,
-        )
-    return 'capped_cylinder', pars
-
-
-@model('cscyl')
-def core_shell_cylinder():
-    pars = dict(
-        scale=1, background=0,
-        core_sld=6, shell_sld=8, solvent_sld=1,
-        radius=45, thickness=25, length=340,
-        theta=30, phi=15,
-        radius_pd=.2, radius_pd_n=1,
-        length_pd=.2, length_pd_n=1,
-        thickness_pd=.2, thickness_pd_n=1,
-        theta_pd=15, theta_pd_n=45,
-        phi_pd=15, phi_pd_n=1,
-        )
-    return 'core_shell_cylinder', pars
-
-
-@model('ell')
-def ellipsoid():
-    pars = dict(
-        scale=1, background=0,
-        sld=6, solvent_sld=1,
-        rpolar=50, requatorial=30,
-        theta=30, phi=15,
-        rpolar_pd=.2, rpolar_pd_n=1,
-        requatorial_pd=.2, requatorial_pd_n=1,
-        theta_pd=15, theta_pd_n=45,
-        phi_pd=15, phi_pd_n=1,
-        )
-    return 'ellipsoid', pars
-
-
-@model('ell3')
-def triaxial_ellipsoid():
-    pars = dict(
-        scale=1, background=0,
-        sld=6, solvent_sld=1,
-        theta=30, phi=15, psi=5,
-        req_minor=25, req_major=36, rpolar=50,
-        req_minor_pd=0, req_minor_pd_n=1,
-        req_major_pd=0, req_major_pd_n=1,
-        rpolar_pd=.2, rpolar_pd_n=1,
-        theta_pd=15, theta_pd_n=45,
-        phi_pd=15, phi_pd_n=1,
-        psi_pd=15, psi_pd_n=1,
-        )
-    return 'triaxial_ellipsoid', pars
-
-@model('sph')
-def sphere():
-    pars = dict(
-        scale=1, background=0,
-        sld=6, solvent_sld=1,
-        radius=120,
-        radius_pd=.2, radius_pd_n=45,
-        )
-    return 'sphere', pars
-
-@model('lam')
-def lamellar():
-    pars = dict(
-        scale=1, background=0,
-        sld=6, solvent_sld=1,
-        thickness=40,
-        thickness_pd= 0.2, thickness_pd_n=40,
-        )
-    return 'lamellar', pars
-
-if __name__ == "__main__":
-    main()
Index: compare.py
===================================================================
--- compare.py	(revision a953943cb9a31b3d6cd9136901b538c74c27179e)
+++ compare.py	(revision 87985ca5989a70a24a7e9ac09a971d544fe16dfe)
@@ -2,10 +2,22 @@
 # -*- coding: utf-8 -*-
 
-from sasmodel import SasModel, load_data, set_beam_stop, plot_data, tic, toc
+import sys
+import math
+
 import numpy as np
 
+from sasmodels.bumps_model import BumpsModel, plot_data, tic
+from sasmodels import gpu, dll
+from sasmodels.convert import revert_model
+
 
 def sasview_model(modelname, **pars):
-    modelname = modelname+"Model"
+    """
+    Load a sasview model given the model name.
+    """
+    # convert model parameters from sasmodel form to sasview form
+    #print "old",sorted(pars.items())
+    modelname, pars = revert_model(modelname, pars)
+    #print "new",sorted(pars.items())
     sans = __import__('sans.models.'+modelname)
     ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None)
@@ -21,187 +33,373 @@
         elif k.endswith("_pd_nsigma"):
             model.dispersion[k[:-10]]['nsigmas'] = v
+        elif k.endswith("_pd_type"):
+            model.dispersion[k[:-8]]['type'] = v
         else:
             model.setParam(k, v)
     return model
 
-
-def sasview_eval(model, data):
-    theory = model.evalDistribution([data.qx_data, data.qy_data])
-    return theory
-
-def plot(data, cpumodel, gpumodel, N, pars):
-
-    model = SasModel(data, gpumodel, dtype='f', **pars)
-    tic()
-    for i in range(N):
-        #pars['scale'] = np.random.rand()
-        model.update()
-        gpu = model.theory()
-    gpu_time = toc()*1000./N
-    print "ocl t=%.1f ms"%gpu_time
-    tic()
-    for i in range(N):
-        cpu = sasview_eval(cpumodel, data)
-    cpu_time = toc()*1000./N
-    print "omp t=%.1f ms"%cpu_time
-
+def load_opencl(modelname, dtype='single'):
+    sasmodels = __import__('sasmodels.models.'+modelname)
+    module = getattr(sasmodels.models, modelname, None)
+    kernel = gpu.load_model(module, dtype=dtype)
+    return kernel
+
+def load_ctypes(modelname, dtype='single'):
+    sasmodels = __import__('sasmodels.models.'+modelname)
+    module = getattr(sasmodels.models, modelname, None)
+    kernel = dll.load_model(module, dtype=dtype)
+    return kernel
+
+def randomize(p, v):
+    """
+    Randomizing parameter.
+
+    Guess the parameter type from name.
+    """
+    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')):
+        return v
+    elif any(s in p for s in ('theta','phi','psi')):
+        # orientation in [-180,180], orientation pd in [0,45]
+        if p.endswith('_pd'):
+            return 45*np.random.rand()
+        else:
+            return 360*np.random.rand() - 180
+    elif 'sld' in p:
+        # sld in in [-0.5,10]
+        return 10.5*np.random.rand() - 0.5
+    elif p.endswith('_pd'):
+        # length pd in [0,1]
+        return np.random.rand()
+    else:
+        # length, scale, background in [0,200]
+        return 200*np.random.rand()
+
+def parlist(pars):
+    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
+
+def suppress_pd(pars):
+    """
+    Suppress theta_pd for now until the normalization is resolved.
+
+    May also suppress complete polydispersity of the model to test
+    models more quickly.
+    """
+    for p in pars:
+        if p.endswith("_pd"): pars[p] = 0
+
+def compare(name, pars, Ncpu, Ngpu, opts, set_pars):
+
+    # Sort out data
+    qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05)
+    if "-1d" in opts:
+        from sasmodels.bumps_model import empty_data1D
+        qmax = math.log10(qmax)
+        data = empty_data1D(np.logspace(qmax-3, qmax, 128))
+        index = slice(None, None)
+    else:
+        from sasmodels.bumps_model import empty_data2D, set_beam_stop
+        data = empty_data2D(np.linspace(-qmax, qmax, 128))
+        set_beam_stop(data, 0.004)
+        index = ~data.mask
+    is2D = hasattr(data, 'qx_data')
+
+    # modelling accuracy is determined by dtype and cutoff
+    dtype = 'double' if '-double' in opts else 'single'
+    cutoff_opts = [s for s in opts if s.startswith('-cutoff')]
+    cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5
+
+    # randomize parameters
+    random_opts = [s for s in opts if s.startswith('-random')]
+    if random_opts:
+        if '=' in random_opts[0]:
+            seed = int(random_opts[0].split('=')[1])
+        else:
+            seed = int(np.random.rand()*10000)
+        np.random.seed(seed)
+        print "Randomize using -random=%i"%seed
+        # Note: the sort guarantees order of calls to random number generator
+        pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items()))
+        # The capped cylinder model has a constraint on its parameters
+        if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
+            pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
+    pars.update(set_pars)
+
+    # parameter selection
+    if '-mono' in opts:
+        suppress_pd(pars)
+    if '-pars' in opts:
+        print "pars",parlist(pars)
+
+    # OpenCl calculation
+    if Ngpu > 0:
+        try:
+            model = load_opencl(name, dtype=dtype)
+        except Exception,exc:
+            print exc
+            print "... trying again with single precision"
+            model = load_opencl(name, dtype='single')
+        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
+        toc = tic()
+        for _ in range(Ngpu):
+            #pars['scale'] = np.random.rand()
+            problem.update()
+            gpu = problem.theory()
+        gpu_time = toc()*1000./Ngpu
+        print "opencl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[index]))
+        #print max(gpu), min(gpu)
+
+    # ctypes/sasview calculation
+    if Ncpu > 0 and "-ctypes" in opts:
+        model = load_ctypes(name, dtype=dtype)
+        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
+        toc = tic()
+        for _ in range(Ncpu):
+            problem.update()
+            cpu = problem.theory()
+        cpu_time = toc()*1000./Ncpu
+        comp = "ctypes"
+        print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
+    elif Ncpu > 0:
+        model = sasview_model(name, **pars)
+        toc = tic()
+        for _ in range(Ncpu):
+            if is2D:
+                cpu = model.evalDistribution([data.qx_data, data.qy_data])
+            else:
+                cpu = model.evalDistribution(data.x)
+        cpu_time = toc()*1000./Ncpu
+        comp = "sasview"
+        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
+
+    # Compare, but only if computing both forms
+    if Ngpu > 0 and Ncpu > 0:
+        #print "speedup %.2g"%(cpu_time/gpu_time)
+        #print "max |gpu/cpu|", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu))
+        #cpu *= max(gpu/cpu)
+        resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu)
+        resid[index] = (gpu - cpu)[index]
+        relerr[index] = resid[index]/cpu[index]
+        print "max(|ocl-%s|)"%comp, max(abs(resid[index]))
+        print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[index]))
+
+    # Plot if requested
+    if '-noplot' in opts: return
     import matplotlib.pyplot as plt
-
-    print "gpu/cpu", max(gpu/cpu)
-    cpu *= max(gpu/cpu)
-    relerr = (gpu - cpu)/cpu
-    print "max(|(ocl-omp)/ocl|)", max(abs(relerr))
-
-
-    plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time)
-    plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time)
-    plt.subplot(133); plot_data(data, relerr); plt.title("relerr x 10^8"); plt.colorbar()
+    if Ncpu > 0:
+        if Ngpu > 0: plt.subplot(131)
+        plot_data(data, cpu, scale='log')
+        plt.title("%s t=%.1f ms"%(comp,cpu_time))
+    if Ngpu > 0:
+        if Ncpu > 0: plt.subplot(132)
+        plot_data(data, gpu, scale='log')
+        plt.title("opencl t=%.1f ms"%gpu_time)
+    if Ncpu > 0 and Ngpu > 0:
+        plt.subplot(133)
+        err = resid if '-abs' in opts else relerr
+        errstr = "abs err" if '-abs' in opts else "rel err"
+        #err,errstr = gpu/cpu,"ratio"
+        plot_data(data, err, scale='linear')
+        plt.title("max %s = %.3g"%(errstr, max(abs(err[index]))))
+        if is2D: plt.colorbar()
     plt.show()
 
-def newlen(N):
-    import sys
-
-    if len(sys.argv) > 1:
-        N = int(sys.argv[1])
-    data = load_data('December/DEC07098.DAT')
-    set_beam_stop(data, 0.004)
-
-    return data, N
-
-def cyl(N=1):
-
-    dtype = "float"
-    pars = dict(
-        scale=.003, radius=64.1, length=66.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1,
-        cyl_theta=90, cyl_phi=0,
-        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
-        length_pd=0.1,length_pd_n=5, length_pd_nsigma=3,
-        cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3,
-        cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3)
-
-    from Models.code_cylinder_f import GpuCylinder as gpumodel
-    model = sasview_model('Cylinder', **pars)
-    data, N = newlen(N)
-
-    plot(data, model, gpumodel, N, pars)
-
-
-def ellipse(N=1):
-
-    pars = dict(scale=.027, radius_a=60, radius_b=180, sldEll=.297e-6, sldSolv=5.773e-6, background=4.9,
-                axis_theta=0, axis_phi=90,
-                radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3,
-                radius_b_pd=0.1, radius_b_pd_n=10, radius_b_pd_nsigma=3,
-                axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3,
-                axis_phi_pd=0.1, axis_phi_pd_n=6, axis_phi_pd_nsigma=3,)
-
-    from Models.code_ellipse import GpuEllipse as gpumodel
-    model = sasview_model('Ellipsoid', **pars)
-
-    data, N = newlen(N)
-    plot(data, model, gpumodel, N, pars)
-
-def coreshell(N=1):
-
-    data, N = newlen(N)
-
-    pars = dict(scale= 1.77881e-06, radius=325, thickness=25, length=34.2709,
-                 core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6,
-                 background=223.827, axis_theta=90, axis_phi=0,
-                 axis_theta_pd=15.8,
-                 radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
-                 length_pd=0.1, length_pd_n=10, length_pd_nsigma=3,
-                 thickness_pd=0.1, thickness_pd_n=5, thickness_pd_nsigma=3,
-                 axis_theta_pd_n=20, axis_theta_pd_nsigma=5,
-                 axis_phi_pd=0.0008748, axis_phi_pd_n=5, axis_phi_pd_nsigma=3,)
-
-    model = sasview_model('CoreShellCylinder', **pars)
-    from Models.code_coreshellcyl_f import GpuCoreShellCylinder as gpumodel
-
-    plot(data, model, gpumodel, N, pars)
-
-def trellipse(N=1):
-
-    data, N = newlen(N)
-
-    pars = dict(scale=0.08, semi_axisA=15, semi_axisB=20, semi_axisC=500,
-                sldEll=7.105e-6, sldSolv=.291e-6,
-                background=5, axis_theta=0, axis_phi=0, axis_psi=0,
-                axis_theta_pd=20, axis_theta_pd_n=10, axis_theta_pd_nsigma=3,
-                axis_phi_pd=.1, axis_phi_pd_n=10, axis_phi_pd_nsigma=3,
-                axis_psi_pd=30, axis_psi_pd_n=5, axis_psi_pd_nsigma=3,
-                semi_axisA_pd=.1, semi_axisA_pd_n=5, semi_axisA_pd_nsigma=3,
-                semi_axisB_pd=.1, semi_axisB_pd_n=5, semi_axisB_pd_nsigma=3,
-                semi_axisC_pd=.1, semi_axisC_pd_n=5, semi_axisC_pd_nsigma=3,)
-
-    model = sasview_model('TriaxialEllipsoid', **pars)
-    from Models.code_triaxialellipse import GpuTriEllipse as gpumodel
-
-    plot(data, model, gpumodel, N, pars)
-
-def lamellar(N=1):
-
-    data, N = newlen(N)
-
-    pars = dict(scale=0.08,
-                bi_thick=19.2946,
-                sld_bi=5.38e-6,sld_sol=7.105e-6,
-                background=0.003,
-                bi_thick_pd= 0.37765, bi_thick_pd_n=40, bi_thick_pd_nsigma=3)
-
-    model = sasview_model('Lamellar', **pars)
-    from Models.code_lamellar import GpuLamellar as gpumodel
-
-    plot(data, model, gpumodel, N, pars)
-
-def cap(N=1):
-
-    data, N = newlen(N)
-
-    pars = dict(scale=.08, rad_cyl=20, rad_cap=40, len_cyl=400,
-                sld_capcyl=1e-6, sld_solv=6.3e-6,
-                background=0, theta=0, phi=0,
-                rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_pd_nsigma=3,
-                rad_cap_pd=.1, rad_cap_pd_n=10, rad_cap_pd_nsigma=3,
-                len_cyl_pd=.1, len_cyl_pd_n=3, len_cyl_pd_nsigma=3,
-                theta_pd=.1, theta_pd_n=3, theta_pd_nsigma=3,
-                phi_pd=.1, phi_pd_n=3, phi_pd_nsigma=3)
-
-
-    model = sasview_model('CappedCylinder', **pars)
-    from Models.code_capcyl import GpuCapCylinder as gpumodel
-
-    plot(data, model, gpumodel, N, pars)
-
+# ===========================================================================
+#
+USAGE="""
+usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val]
+
+Compare the speed and value for a model between the SasView original and the
+OpenCL rewrite.
+
+model is the name of the model to compare (see below).
+Nopencl is the number of times to run the OpenCL model (default=5)
+Nsasview is the number of times to run the Sasview model (default=1)
+
+Options (* for default):
+
+    -plot*/-noplot plots or suppress the plot of the model
+    -single/-double uses double precision for comparison
+    -lowq/-midq/-highq use q values up to 0.05, 0.2 or 1.0
+    -2d/-1d uses 1d or 2d random data
+    -preset/-random[=seed] preset or random parameters
+    -poly/-mono force monodisperse/polydisperse
+    -sasview/-ctypes whether cpu is tested using sasview or ctypes
+    -cutoff=1e-5/value cutoff for including a point in polydispersity
+    -nopars/-pars prints the parameter set or not
+    -rel/-abs plot relative or absolute error
+
+Key=value pairs allow you to set specific values to any of the model
+parameters.
+
+Available models:
+
+    %s
+"""
+
+VALID_OPTIONS = [
+    'plot','noplot',
+    'single','double',
+    'lowq','midq','highq',
+    '2d','1d',
+    'preset','random',
+    'poly','mono',
+    'sasview','ctypes',
+    'nopars','pars',
+    'rel','abs',
+    ]
+
+def main():
+    opts = [arg for arg in sys.argv[1:] if arg.startswith('-')]
+    args = [arg for arg in sys.argv[1:] if not arg.startswith('-')]
+    models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' '))
+                           for k,v in sorted(MODELS.items()))
+    if len(args) == 0:
+        print(USAGE%models)
+        sys.exit(1)
+    if args[0] not in MODELS:
+        print "Model %r not available. Use one of:\n    %s"%(args[0],models)
+        sys.exit(1)
+
+    valid_opts = set(VALID_OPTIONS)
+    invalid = [o[1:] for o in opts
+               if o[1:] not in valid_opts
+                  and not o.startswith('-cutoff=')
+                  and not o.startswith('-random=')]
+    if invalid:
+        print "Invalid options: %s"%(", ".join(invalid))
+        sys.exit(1)
+
+    name, pars = MODELS[args[0]]()
+    Nopencl = int(args[1]) if len(args) > 1 else 5
+    Nsasview = int(args[2]) if len(args) > 3 else 1
+
+    # Fill in default polydispersity parameters
+    pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd'))
+    for p in pds:
+        if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3
+        if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian"
+
+    # Fill in parameters given on the command line
+    set_pars = {}
+    for arg in args[3:]:
+        k,v = arg.split('=')
+        if k not in pars:
+            # extract base name without distribution
+            s = set(p.split('_pd')[0] for p in pars)
+            print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s)))
+            sys.exit(1)
+        set_pars[k] = float(v) if not v.endswith('type') else v
+
+    compare(name, pars, Nsasview, Nopencl, opts, set_pars)
+
+# ===========================================================================
+#
+
+MODELS = {}
+def model(name):
+    def gather_function(fn):
+        MODELS[name] = fn
+        return fn
+    return gather_function
+
+
+@model('cyl')
+def cylinder():
+    pars = dict(
+        scale=1, background=0,
+        sld=6, solvent_sld=1,
+        #radius=5, length=20,
+        radius=260, length=290,
+        theta=30, phi=0,
+        radius_pd=.2, radius_pd_n=1,
+        length_pd=.2,length_pd_n=1,
+        theta_pd=15, theta_pd_n=45,
+        phi_pd=15, phi_pd_n=1,
+        )
+    return 'cylinder', pars
+
+@model('capcyl')
+def capped_cylinder():
+    pars = dict(
+        scale=1, background=0,
+        sld=6, solvent_sld=1,
+        radius=260, cap_radius=290, length=290,
+        theta=30, phi=15,
+        radius_pd=.2, radius_pd_n=1,
+        cap_radius_pd=.2, cap_radius_pd_n=1,
+        length_pd=.2, length_pd_n=1,
+        theta_pd=15, theta_pd_n=45,
+        phi_pd=15, phi_pd_n=1,
+        )
+    return 'capped_cylinder', pars
+
+
+@model('cscyl')
+def core_shell_cylinder():
+    pars = dict(
+        scale=1, background=0,
+        core_sld=6, shell_sld=8, solvent_sld=1,
+        radius=45, thickness=25, length=340,
+        theta=30, phi=15,
+        radius_pd=.2, radius_pd_n=1,
+        length_pd=.2, length_pd_n=1,
+        thickness_pd=.2, thickness_pd_n=1,
+        theta_pd=15, theta_pd_n=45,
+        phi_pd=15, phi_pd_n=1,
+        )
+    return 'core_shell_cylinder', pars
+
+
+@model('ell')
+def ellipsoid():
+    pars = dict(
+        scale=1, background=0,
+        sld=6, solvent_sld=1,
+        rpolar=50, requatorial=30,
+        theta=30, phi=15,
+        rpolar_pd=.2, rpolar_pd_n=1,
+        requatorial_pd=.2, requatorial_pd_n=1,
+        theta_pd=15, theta_pd_n=45,
+        phi_pd=15, phi_pd_n=1,
+        )
+    return 'ellipsoid', pars
+
+
+@model('ell3')
+def triaxial_ellipsoid():
+    pars = dict(
+        scale=1, background=0,
+        sld=6, solvent_sld=1,
+        theta=30, phi=15, psi=5,
+        req_minor=25, req_major=36, rpolar=50,
+        req_minor_pd=0, req_minor_pd_n=1,
+        req_major_pd=0, req_major_pd_n=1,
+        rpolar_pd=.2, rpolar_pd_n=1,
+        theta_pd=15, theta_pd_n=45,
+        phi_pd=15, phi_pd_n=1,
+        psi_pd=15, psi_pd_n=1,
+        )
+    return 'triaxial_ellipsoid', pars
+
+@model('sph')
+def sphere():
+    pars = dict(
+        scale=1, background=0,
+        sld=6, solvent_sld=1,
+        radius=120,
+        radius_pd=.2, radius_pd_n=45,
+        )
+    return 'sphere', pars
+
+@model('lam')
+def lamellar():
+    pars = dict(
+        scale=1, background=0,
+        sld=6, solvent_sld=1,
+        thickness=40,
+        thickness_pd= 0.2, thickness_pd_n=40,
+        )
+    return 'lamellar', pars
 
 if __name__ == "__main__":
-    coreshell(1)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+    main()
Index: fit.py
===================================================================
--- fit.py	(revision 4001d6ef0c7f2735a7194f740e6ce959bca030e5)
+++ fit.py	(revision 87985ca5989a70a24a7e9ac09a971d544fe16dfe)
@@ -2,90 +2,80 @@
 # -*- coding: utf-8 -*-
 
+import sys
 from bumps.names import *
-from sasmodel import SasModel, load_data, set_beam_stop, set_half, set_top
-from Models.code_capcyl import GpuCapCylinder
-from Models.code_coreshellcyl_f import GpuCoreShellCylinder
-from Models.code_cylinder_f import GpuCylinder
-#from Models.code_cylinder import GpuCylinder, OneDGpuCylinder
-from Models.code_ellipse_f import GpuEllipse
-from Models.code_lamellar import GpuLamellar
-from Models.code_triaxialellipse import GpuTriEllipse
+from sasmodels import bumps_model as sas
 
 """ IMPORT THE DATA USED """
-#data = load_data('December/Tangential/Sector0/DEC07133.ABS')
-data = load_data('December/DEC07098.DAT')
+radial_data = sas.load_data('December/DEC07267.DAT')
+sas.set_beam_stop(radial_data, 0.00669, outer=0.025)
+sas.set_top(radial_data, -.0185)
 
-""" 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_half(data, 'right')
+tan_data = sas.load_data('December/DEC07266.DAT')
+sas.set_beam_stop(tan_data, 0.00669, outer=0.025)
+sas.set_top(tan_data, -.0185)
+#sas.set_half(tan_data, 'right')
 
+name = "ellipsoid" if len(sys.argv) < 2 else sys.argv[1]
+section = "radial" if len(sys.argv) < 3 else sys.argv[2]
+data = radial_data if section is not "tangent" else tan_data
+kernel = sas.load_model(name, dtype="single")
 
-if 0:
-    model = SasModel(data, OneDGpuCylinder,
-    scale=0.0013, radius=105, length=1000,
-    background=21, sldCyl=.291e-6,sldSolv=7.105e-6,
-    radius_pd=0.1,radius_pd_n=10,radius_pd_nsigma=0,
-    length_pd=0.1,length_pd_n=5,length_pd_nsigma=0,
-    bolim=0.0, uplim=90) #bottom limit, upper limit of angle integral
-
-
-if 0:
-    model = SasModel(data, GpuEllipse,
+if name == "ellipsoid":
+    model = sas.BumpsModel(data, kernel,
     scale=0.08,
-    radius_a=15, radius_b=800,
-    sldEll=.291e-6, sldSolv=7.105e-6,
+    rpolar=15, requatorial=800,
+    sld=.291, solvent_sld=7.105,
     background=0,
-    axis_theta=90, axis_phi=0,
-    axis_theta_pd=15, axis_theta_pd_n=40, axis_theta_pd_nsigma=3,
-    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=3,
-    dtype='float32')
+    theta=90, phi=0,
+    theta_pd=15, theta_pd_n=40, theta_pd_nsigma=3,
+    rpolar_pd=0.222296, rpolar_pd_n=1, rpolar_pd_nsigma=0,
+    requatorial_pd=.000128, requatorial_pd_n=1, requatorial_pd_nsigma=0,
+    phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3,
+    )
 
 
     # SET THE FITTING PARAMETERS
-    #model.radius_a.range(15, 1000)
-    #model.radius_b.range(15, 1000)
-    #model.axis_theta_pd.range(0, 360)
+    #model.rpolar.range(15, 1000)
+    #model.requatorial.range(15, 1000)
+    #model.theta_pd.range(0, 360)
     #model.background.range(0,1000)
     model.scale.range(0, 1)
 
 
-if 0:
-    model = SasModel(data, GpuLamellar,
+elif name == "lamellar":
+    model = sas.BumpsModel(data, kernel,
     scale=0.08,
-    bi_thick=19.2946,
-    sld_bi=5.38e-6,sld_sol=7.105e-6,
+    thickness=19.2946,
+    sld=5.38,sld_sol=7.105,
     background=0.003,
-    bi_thick_pd= 0.37765, bi_thick_pd_n=10, bi_thick_pd_nsigma=3,
-    dtype='float32')
+    thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3,
+    )
 
     # SET THE FITTING PARAMETERS
-    #model.bi_thick.range(0, 1000)
+    #model.thickness.range(0, 1000)
     #model.scale.range(0, 1)
-    #model.bi_thick_pd.range(0, 1000)
+    #model.thickness_pd.range(0, 1000)
     #model.background.range(0, 1000)
-    model.sld_bi.range(0, 1)
+    model.sld.range(0, 1)
 
 
-if 0:
+elif name == "cylinder":
     """
     pars = dict(scale=0.0023, radius=92.5, length=798.3,
-        sldCyl=.29e-6, sldSolv=7.105e-6, background=5,
-        cyl_theta=0, cyl_phi=0,
-        cyl_theta_pd=22.11, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3,
+        sld=.29, solvent_sld=7.105, background=5,
+        theta=0, phi=0,
+        theta_pd=22.11, theta_pd_n=5, theta_pd_nsigma=3,
         radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3,
         length_pd=0.493, length_pd_n=10, length_pd_nsigma=3,
-        cyl_phi_pd=0, cyl_phi_pd_n=5, cyl_phi_pd_nsigma=3,)
+        phi_pd=0, phi_pd_n=5, phi_pd_nsigma=3,)
         """
     pars = dict(
-        scale=.01, radius=64.1, length=66.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1,
-        cyl_theta=90, cyl_phi=0,
+        scale=.01, radius=64.1, length=66.96, sld=.291, solvent_sld=5.77, background=.1,
+        theta=90, phi=0,
         radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
         length_pd=0.1,length_pd_n=5, length_pd_nsigma=3,
-        cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3,
-        cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3)
-    model = SasModel(data, GpuCylinder, dtype="float32", **pars)
+        theta_pd=0.1, theta_pd_n=5, theta_pd_nsigma=3,
+        phi_pd=0.1, phi_pd_n=10, phi_pd_nsigma=3)
+    model = sas.BumpsModel(data, kernel, **pars)
 
 
@@ -94,26 +84,26 @@
     model.radius.range(1, 500)
     model.length.range(1, 5000)
-    #model.cyl_theta.range(-90,100)
-    #model.cyl_theta_pd.range(0, 90)
-    #model.cyl_theta_pd_n = model.cyl_theta_pd + 5
+    #model.theta.range(-90,100)
+    #model.theta_pd.range(0, 90)
+    #model.theta_pd_n = model.theta_pd + 5
     #model.radius_pd.range(0, 90)
     #model.length_pd.range(0, 90)
     model.scale.range(0, 1)
     #model.background.range(0, 100)
-    #model.sldCyl.range(0, 1)
+    #model.sld.range(0, 1)
 
 
-if 1:
-    model = SasModel(data, GpuCoreShellCylinder,
+elif name == "core_shell_cylinder":
+    model = sas.BumpsModel(data, kernel,
     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, axis_theta=0, axis_phi=0,
+    core_sld=7.105, shell_sld=.291, solvent_sld=7.105,
+    background=0, theta=0, phi=0,
 
     radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3,
     length_pd=0.26, length_pd_n=10, length_pd_nsigma=3,
     thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1,
-    axis_theta_pd=1, axis_theta_pd_n=10, axis_theta_pd_nsigma=3,
-    axis_phi_pd=0.1, axis_phi_pd_n=1, axis_phi_pd_nsigma=1,
-    dtype='float32')
+    theta_pd=1, theta_pd_n=10, theta_pd_nsigma=3,
+    phi_pd=0.1, phi_pd_n=1, phi_pd_nsigma=1,
+    )
 
     # SET THE FITTING PARAMETERS
@@ -122,8 +112,8 @@
     #model.thickness.range(18, 38)
     #model.thickness_pd.range(0, 1)
-    #model.axis_phi.range(0, 90)
+    #model.phi.range(0, 90)
     #model.radius_pd.range(0, 1)
     #model.length_pd.range(0, 1)
-    #model.axis_theta_pd.range(0, 360)
+    #model.theta_pd.range(0, 360)
     #model.background.range(0,5)
     model.scale.range(0, 1)
@@ -131,35 +121,35 @@
 
 
-if 0:
-    model = SasModel(data, GpuCapCylinder,
-    scale=.08, rad_cyl=20, rad_cap=40, len_cyl=400,
-    sld_capcyl=1e-6, sld_solv=6.3e-6,
+elif name == "capped_cylinder":
+    model = sas.BumpsModel(data, kernel,
+    scale=.08, radius=20, cap_radius=40, length=400,
+    sld_capcyl=1, sld_solv=6.3,
     background=0, theta=0, phi=0,
-    rad_cyl_pd=.1, rad_cyl_pd_n=5, rad_cyl_pd_nsigma=3,
-    rad_cap_pd=.1, rad_cap_pd_n=5, rad_cap_pd_nsigma=3,
-    len_cyl_pd=.1, len_cyl_pd_n=1, len_cyl_pd_nsigma=0,
+    radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3,
+    cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3,
+    length_pd=.1, length_pd_n=1, length_pd_nsigma=0,
     theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0,
     phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0,
-    dtype='float32')
+    dtype=dtype)
 
     model.scale.range(0, 1)
 
 
-if 0:
-    model = SasModel(data, GpuTriEllipse,
-    scale=0.08, axisA=15, axisB=20, axisC=500,
-    sldEll=7.105e-6, sldSolv=.291e-6,
+elif name == "triaxial_ellipsoid":
+    model = sas.BumpsModel(data, kernel,
+    scale=0.08, req_minor=15, req_major=20, rpolar=500,
+    sldEll=7.105, solvent_sld=.291,
     background=5, theta=0, phi=0, psi=0,
     theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3,
     phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0,
     psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0,
-    axisA_pd=.1, axisA_pd_n=1, axisA_pd_nsigma=0,
-    axisB_pd=.1, axisB_pd_n=1, axisB_pd_nsigma=0,
-    axisC_pd=.1, axisC_pd_n=1, axisC_pd_nsigma=0, dtype='float32')
+    req_minor_pd=.1, req_minor_pd_n=1, req_minor_pd_nsigma=0,
+    req_major_pd=.1, req_major_pd_n=1, req_major_pd_nsigma=0,
+    rpolar_pd=.1, rpolar_pd_n=1, rpolar_pd_nsigma=0, dtype=dtype)
 
     # SET THE FITTING PARAMETERS
-    model.axisA.range(15, 1000)
-    model.axisB.range(15, 1000)
-    #model.axisC.range(15, 1000)
+    model.req_minor.range(15, 1000)
+    model.req_major.range(15, 1000)
+    #model.rpolar.range(15, 1000)
     #model.background.range(0,1000)
     #model.theta_pd.range(0, 360)
@@ -167,7 +157,13 @@
     #model.psi_pd.range(0, 360)
 
+else:
+    print "No parameters for %s"%name
+    sys.exit(1)
 
+if section is "both":
+   tan_model = sas.BumpsModel(tan_data, model.model, model.parameters())
+   tan_model.phi = model.phi - 90
+   problem = FitProblem([model, tan_model])
+else:
+   problem = FitProblem(model)
 
-
-problem = FitProblem(model)
-
Index: t2.py
===================================================================
--- fit2.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ 	(revision )
@@ -1,53 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-from bumps.names import *
-from Models.code_coreshellcyl import GpuCoreShellCylinder
-from sasmodel import SasModel, load_data, set_beam_stop, set_top
-
-
-radial_data = load_data('December/DEC07267.DAT')
-set_beam_stop(radial_data, 0.00669, outer=0.025)
-set_top(radial_data, -.0185)
-
-tan_data = load_data('December/DEC07266.DAT')
-set_beam_stop(tan_data, 0.00669, outer=0.025)
-set_top(tan_data, -.0185)
-
-
-
-
-dtype='float32'
-radial = SasModel(radial_data,
-                 GpuCoreShellCylinder,
-                 scale= 3.75e-7, radius=378, thickness=30, length=1806,
-                 core_sld=7.105e-6, shell_sld=.291e-6, solvent_sld=7.105e-6,
-                 background=0.2, axis_theta=0, axis_phi=90,
-
-                 radius_pd=0.26, radius_pd_n=20, radius_pd_nsigma=3,
-                 length_pd=0.26, length_pd_n=20, length_pd_nsigma=3,
-                 thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=0,
-                 axis_theta_pd=1, axis_theta_pd_n=10, axis_theta_pd_nsigma=3,
-                 axis_phi_pd=0.1, axis_phi_pd_n=1, axis_phi_pd_nsigma=0,
-                 dtype='float')
-tan = SasModel(tan_data,
-                  GpuCoreShellCylinder, dtype=dtype,
-                  **radial.parameters())
-
-radial.radius.range(15, 1000)
-radial.length.range(0, 2500)
-#radial.thickness.range(18, 38)
-#radial.thickness_pd.range(0, 1)
-#radial.axis_phi.range(0, 90)
-#radial.radius_pd.range(0, 1)
-#radial.length_pd.range(0, 1)
-#radial.axis_theta_pd.range(0, 360)
-#radial.background.range(0,5)
-#radial.scale.range(0, 1)
-radial.axis_phi = tan.axis_phi + 90
-
-
-problem = FitProblem([radial,tan])
-
-
-
Index: lti-fit.py
===================================================================
--- multi-fit.py	(revision a42fec0a6299c921e6b9422a1e4fa8aa60d78022)
+++ 	(revision )
@@ -1,61 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-from bumps.names import *
-from Models.code_lamellar import GpuLamellar
-from Models.code_ellipse import GpuEllipse
-from Models.code_cylinder import GpuCylinder
-from multisasmodels import SasModel, load_data, set_beam_stop, set_half
-import numpy as np
-
-data = load_data('December/DEC07238.DAT')
-set_beam_stop(data, 0.0052)#, outer=0.025)
-#set_half(data, 'left')
-
-truth = SasModel(data, GpuCylinder,
-scale=0.08, radius=46, length=719,
-sldCyl=.291e-6, sldSolv=5.77e-6, background=0,
-cyl_theta=90, cyl_phi=0,
-cyl_theta_pd=23, cyl_theta_pd_n=40, cyl_theta_pd_nsigma=3,
-radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
-length_pd=0.1, length_pd_n=10, length_pd_nsigma=3,
-cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3,
-dtype='float')
-
-
-#model.radius_a.range(15, 1000)
-#model.radius_b.range(15, 1000)
-#model.axis_theta_pd.range(0, 360)
-#model.background.range(0,1000)
-#truth.scale.range(0, 1)
-
-arrayone = truth.theory()
-
-
-lies = SasModel(data, GpuCylinder,
-scale=0.08, radius=46, length=719,
-sldCyl=.291e-6, sldSolv=5.77e-6, background=0,
-cyl_theta=90, cyl_phi=0,
-cyl_theta_pd=23, cyl_theta_pd_n=40, cyl_theta_pd_nsigma=3,
-radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
-length_pd=0.1, length_pd_n=10, length_pd_nsigma=3,
-cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3,
-dtype='float')
-
-
-
-#modeltwo.bi_thick.range(0, 1000)
-#modeltwo.scale.range(0, 1)
-#model.bi_thick_pd.range(0, 1000)
-#model.background.range(0, 1000)
-
-
-arraytwo = lies.theory()
-
-
-a = np.add(np.multiply(arrayone, .08), np.multiply(arraytwo, .08))
-truth.set_result(a)
-
-
-problem = FitProblem(truth)
-
Index: ltisasmodels.py
===================================================================
--- multisasmodels.py	(revision d772f5d5cfc1cb59c225465b709386ac13da531b)
+++ 	(revision )
@@ -1,157 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import sys
-import numpy as np
-import pyopencl as cl
-from bumps.names import Parameter
-from sans.dataloader.loader import Loader
-from sans.dataloader.manipulations import Ringcut, Boxcut
-
-
-def load_data(filename):
-    loader = Loader()
-    data = loader.load(filename)
-    if data is None:
-        raise IOError("Data %r could not be loaded"%filename)
-    return data
-
-
-def set_beam_stop(data, radius, outer=None):
-    data.mask = Ringcut(0, radius)(data)
-    if outer is not None:
-        data.mask += Ringcut(outer,np.inf)(data)
-
-def set_half(data, half):
-    if half == 'left':
-        data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
-    if half == 'right':
-        data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
-
-
-
-def plot_data(data, iq, vmin=None, vmax=None):
-    from numpy.ma import masked_array
-    import matplotlib.pyplot as plt
-    img = masked_array(iq, data.mask)
-    xmin, xmax = min(data.qx_data), max(data.qx_data)
-    ymin, ymax = min(data.qy_data), max(data.qy_data)
-    plt.imshow(img.reshape(128,128),
-               interpolation='nearest', aspect=1, origin='upper',
-               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
-
-def plot_result(data, theory, view='linear'):
-    import matplotlib.pyplot as plt
-    from numpy.ma import masked_array, masked
-    plt.subplot(1, 3, 1)
-    #print "not a number",sum(np.isnan(data.data))
-    #data.data[data.data<0.05] = 0.5
-    mdata = masked_array(data.data, data.mask)
-    mdata[np.isnan(mdata)] = masked
-    if view is 'log':
-        mdata[mdata <= 0] = masked
-        mdata = np.log10(mdata)
-        mtheory = masked_array(np.log10(theory), mdata.mask)
-    else:
-        mtheory = masked_array(theory, mdata.mask)
-    mresid = masked_array((theory-data.data)/data.err_data, data.mask)
-    vmin = min(mdata.min(), mtheory.min())
-    vmax = max(mdata.max(), mtheory.max())
-
-    plot_data(data, mdata, vmin=vmin, vmax=vmax)
-    plt.colorbar()
-    plt.subplot(1, 3, 2)
-    plot_data(data, mtheory, vmin=vmin, vmax=vmax)
-    plt.colorbar()
-    plt.subplot(1, 3, 3)
-    plot_data(data, mresid)
-    plt.colorbar()
-
-
-def demo():
-    data = load_data('JUN03289.DAT')
-    set_beam_stop(data, 0.004)
-    plot_data(data)
-    import matplotlib.pyplot as plt; plt.show()
-
-
-GPU_CONTEXT = None
-GPU_QUEUE = None
-def card():
-    global GPU_CONTEXT, GPU_QUEUE
-    if GPU_CONTEXT is None:
-        GPU_CONTEXT = cl.create_some_context()
-        GPU_QUEUE = cl.CommandQueue(GPU_CONTEXT)
-    return GPU_CONTEXT, GPU_QUEUE
-
-
-class SasModel(object):
-    def __init__(self, data, model, dtype='float32', **kw):
-        self.__dict__['_parameters'] = {}
-        self.index = (data.mask==0) & (~np.isnan(data.data))
-        self.iq = data.data[self.index]
-        self.diq = data.err_data[self.index]
-        self.data = data
-        self.qx = data.qx_data
-        self.qy = data.qy_data
-        self.gpu = model(self.qx, self.qy, dtype=dtype)
-        pd_pars = set(base+attr for base in model.PD_PARS for attr in ('_pd','_pd_n','_pd_nsigma'))
-        total_pars = set(model.PARS.keys()) | pd_pars
-        extra_pars = set(kw.keys()) - total_pars
-        if extra_pars:
-            raise TypeError("unexpected parameters %s"%(str(extra_pars,)))
-        pars = model.PARS.copy()
-        pars.update((base+'_pd', 0) for base in model.PD_PARS)
-        pars.update((base+'_pd_n', 35) for base in model.PD_PARS)
-        pars.update((base+'_pd_nsigma', 3) for base in model.PD_PARS)
-        pars.update(kw)
-        self._parameters = dict((k, Parameter.default(v, name=k)) for k, v in pars.items())
-
-    def set_result(self, result):
-        self.result = result
-        return self.result
-
-    def get_result(self):
-        return self.result
-
-    def numpoints(self):
-        return len(self.iq)
-
-    def parameters(self):
-        return self._parameters
-
-    def __getattr__(self, par):
-        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())
-        result = self.gpu.eval(pars)
-        return result
-
-    def residuals(self):
-        #if np.any(self.err ==0): print "zeros in err"
-        return (self.get_result()[self.index]-self.iq)/self.diq
-
-    def nllf(self):
-        R = self.residuals()
-        #if np.any(np.isnan(R)): print "NaN in residuals"
-        return 0.5*np.sum(R**2)
-
-    def __call__(self):
-        return 2*self.nllf()/self.dof
-
-    def plot(self, view='log'):
-        plot_result(self.data, self.get_result(), view=view)
-
-    def save(self, basename):
-        pass
-
-    def update(self):
-        pass
-
Index: smodel.py
===================================================================
--- sasmodel.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
+++ 	(revision )
@@ -1,232 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import datetime
-import numpy as np
-
-
-TIC = None
-def tic():
-    global TIC
-    TIC = datetime.datetime.now()
-
-def toc():
-    now = datetime.datetime.now()
-    return (now-TIC).total_seconds()
-
-def load_data(filename):
-    from sans.dataloader.loader import Loader
-    loader = Loader()
-    data = loader.load(filename)
-    if data is None:
-        raise IOError("Data %r could not be loaded"%filename)
-    return data
-
-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(x) (x##f)
-#define real float
-"""
-    else:
-        header = """\
-#pragma OPENCL EXTENSION cl_khr_fp64: enable
-#define REAL(x) (x)
-#define real double
-"""
-    return header+src, qx, qy
-
-def set_precision_1d(src, q, dtype):
-    q = np.ascontiguousarray(q, 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, q
-
-def set_beam_stop(data, radius, outer=None):
-    from sans.dataloader.manipulations import Ringcut
-    if hasattr(data, 'qx_data'):
-        data.mask = Ringcut(0, radius)(data)
-        if outer is not None:
-            data.mask += Ringcut(outer,np.inf)(data)
-    else:
-        data.mask = (data.x>=radius)
-        if outer is not None:
-            data.mask &= (data.x<outer)
-
-def set_half(data, half):
-    from sans.dataloader.manipulations import Boxcut
-    if half == 'right':
-        data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
-    if half == 'left':
-        data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
-
-def set_top(data, max):
-    from sans.dataloader.manipulations import Boxcut
-    data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data)
-
-def plot_data(data, iq, vmin=None, vmax=None):
-    from numpy.ma import masked_array
-    import matplotlib.pyplot as plt
-    img = masked_array(iq, data.mask)
-    xmin, xmax = min(data.qx_data), max(data.qx_data)
-    ymin, ymax = min(data.qy_data), max(data.qy_data)
-    plt.imshow(img.reshape(128,128),
-               interpolation='nearest', aspect=1, origin='upper',
-               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
-
-def plot_result2D(data, theory, view='linear'):
-    import matplotlib.pyplot as plt
-    from numpy.ma import masked_array, masked
-    #print "not a number",sum(np.isnan(data.data))
-    #data.data[data.data<0.05] = 0.5
-    mdata = masked_array(data.data, data.mask)
-    mdata[np.isnan(mdata)] = masked
-    if view is 'log':
-        mdata[mdata <= 0] = masked
-        mdata = np.log10(mdata)
-        mtheory = masked_array(np.log10(theory), mdata.mask)
-    else:
-        mtheory = masked_array(theory, mdata.mask)
-    mresid = masked_array((theory-data.data)/data.err_data, data.mask)
-    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)
-    plot_data(data, mdata, vmin=vmin, vmax=vmax)
-    plt.colorbar()
-    plt.subplot(1, 3, 2)
-    plot_data(data, mtheory, vmin=vmin, vmax=vmax)
-    plt.colorbar()
-    plt.subplot(1, 3, 3)
-    print abs(mresid).max()
-    plot_data(data, mresid)
-    plt.colorbar()
-
-
-def plot_result1D(data, theory, view='linear'):
-    import matplotlib.pyplot as plt
-    from numpy.ma import masked_array, masked
-    #print "not a number",sum(np.isnan(data.y))
-    #data.y[data.y<0.05] = 0.5
-    mdata = masked_array(data.y, data.mask)
-    mdata[np.isnan(mdata)] = masked
-    if view is 'log':
-        mdata[mdata <= 0] = masked
-    mtheory = masked_array(theory, mdata.mask)
-    mresid = masked_array((theory-data.y)/data.dy, mdata.mask)
-
-    plt.subplot(1,2,1)
-    plt.errorbar(data.x, mdata, yerr=data.dy)
-    plt.plot(data.x, mtheory, '-', hold=True)
-    plt.yscale(view)
-    plt.subplot(1, 2, 2)
-    plt.plot(data.x, mresid, 'x')
-    #plt.axhline(1, color='black', ls='--',lw=1, hold=True)
-    #plt.axhline(0, color='black', lw=1, hold=True)
-    #plt.axhline(-1, color='black', ls='--',lw=1, hold=True)
-
-
-
-GPU_CONTEXT = None
-GPU_QUEUE = None
-def card():
-    global GPU_CONTEXT, GPU_QUEUE
-    if GPU_CONTEXT is None:
-        import pyopencl as cl
-        GPU_CONTEXT = cl.create_some_context()
-        GPU_QUEUE = cl.CommandQueue(GPU_CONTEXT)
-    return GPU_CONTEXT, GPU_QUEUE
-
-
-class SasModel(object):
-    def __init__(self, data, model, dtype='float32', **kw):
-        from bumps.names import Parameter
-
-        self.__dict__['_parameters'] = {}
-        #self.name = data.filename
-        self.is2D = hasattr(data,'qx_data')
-        self.data = data
-        if self.is2D:
-            self.index = (data.mask==0) & (~np.isnan(data.data))
-            self.iq = data.data[self.index]
-            self.diq = data.err_data[self.index]
-            self.qx = data.qx_data
-            self.qy = data.qy_data
-            self.gpu = model(self.qx, self.qy, dtype=dtype)
-        else:
-            self.index = (data.mask==0) & (~np.isnan(data.y))
-            self.iq = data.y[self.index]
-            self.diq = data.dy[self.index]
-            self.q = data.x
-            self.gpu = model(self.q, dtype=dtype)
-        pd_pars = set(base+attr for base in model.PD_PARS for attr in ('_pd','_pd_n','_pd_nsigma'))
-        total_pars = set(model.PARS.keys()) | pd_pars
-        extra_pars = set(kw.keys()) - total_pars
-        if extra_pars:
-            raise TypeError("unexpected parameters %s"%(str(extra_pars,)))
-        pars = model.PARS.copy()
-        pars.update((base+'_pd', 0) for base in model.PD_PARS)
-        pars.update((base+'_pd_n', 35) for base in model.PD_PARS)
-        pars.update((base+'_pd_nsigma', 3) for base in model.PD_PARS)
-        pars.update(kw)
-        for k,v in pars.items():
-            setattr(self, k, Parameter.default(v, name=k))
-        self._parameter_names = set(pars.keys())
-        self.update()
-
-    def update(self):
-        self._cache = {}
-
-    def numpoints(self):
-        return len(self.iq)
-
-    def parameters(self):
-        return dict((k,getattr(self,k)) for k in self._parameter_names)
-
-    def theory(self):
-        if 'theory' not in self._cache:
-            pars = dict((k,getattr(self,k).value) for k in self._parameter_names)
-            #print pars
-            self._cache['theory'] = self.gpu.eval(pars)
-        return self._cache['theory']
-
-    def residuals(self):
-        #if np.any(self.err ==0): print "zeros in err"
-        return (self.theory()[self.index]-self.iq)/self.diq
-
-    def nllf(self):
-        R = self.residuals()
-        #if np.any(np.isnan(R)): print "NaN in residuals"
-        return 0.5*np.sum(R**2)
-
-    def __call__(self):
-        return 2*self.nllf()/self.dof
-
-    def plot(self, view='log'):
-        if self.is2D:
-            plot_result2D(self.data, self.theory(), view=view)
-        else:
-            plot_result1D(self.data, self.theory(), view=view)
-
-    def save(self, basename):
-        pass
-
-def demo():
-    data = load_data('DEC07086.DAT')
-    set_beam_stop(data, 0.004)
-    plot_data(data)
-    import matplotlib.pyplot as plt; plt.show()
-
-if __name__ == "__main__":
-    demo()
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/bumps_model.py	(revision 87985ca5989a70a24a7e9ac09a971d544fe16dfe)
@@ -6,4 +6,20 @@
 
 import numpy as np
+
+try:
+    from .gpu import load_model as _loader
+except ImportError,exc:
+    warnings.warn(str(exc))
+    warnings.warn("OpenCL not available --- using ctypes instead")
+    from .dll import load_model as _loader
+
+def load_model(modelname, dtype='single'):
+    """
+    Load model by name.
+    """
+    sasmodels = __import__('sasmodels.models.'+modelname)
+    module = getattr(sasmodels.models, modelname, None)
+    model = _loader(module, dtype=dtype)
+    return model
 
 
@@ -235,6 +251,9 @@
         partype = model.info['partype']
 
+        # remember inputs so we can inspect from outside
+        self.data = data
+        self.model = model
+
         # interpret data
-        self.data = data
         if hasattr(data, 'qx_data'):
             self.index = (data.mask==0) & (~np.isnan(data.data))
Index: smodels/models/cylinder_onefile.py
===================================================================
--- sasmodels/models/cylinder_onefile.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ 	(revision )
@@ -1,211 +1,0 @@
-# This is how a model would look if we were to stick it all into one file.
-r"""
-CylinderModel
-=============
-
-This model provides the form factor for a right circular cylinder with uniform
-scattering length density. The form factor is normalized by the particle volume.
-
-For information about polarised and magnetic scattering, click here_.
-
-Definition
-----------
-
-The output of the 2D scattering intensity function for oriented cylinders is
-given by (Guinier, 1955)
-
-.. math::
-
-    P(q,\alpha) = \frac{\text{scale}}{V}f^2(q) + \text{bkg}
-
-where
-
-.. math::
-
-    f(q) = 2 (\Delta \rho) V
-           \frac{\sin (q L/2 \cos \alpha)}{q L/2 \cos \alpha}
-           \frac{J_1 (q r \sin \alpha)}{q r \sin \alpha}
-
-and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V$
-is the volume of the cylinder, $L$ is the length of the cylinder, $r$ is the
-radius of the cylinder, and $d\rho$ (contrast) is the scattering length density
-difference between the scatterer and the solvent. $J_1$ is the first order
-Bessel function.
-
-To provide easy access to the orientation of the cylinder, we define the
-axis of the cylinder using two angles $\theta$ and $\phi$. Those angles
-are defined in Figure :num:`figure #CylinderModel-orientation`.
-
-.. _CylinderModel-orientation:
-
-.. figure:: img/image061.JPG
-
-    Definition of the angles for oriented cylinders.
-
-.. figure:: img/image062.JPG
-
-    Examples of the angles for oriented pp against the detector plane.
-
-NB: The 2nd virial coefficient of the cylinder is calculated based on the
-radius and length values, and used as the effective radius for $S(Q)$
-when $P(Q) \cdot S(Q)$ is applied.
-
-The returned value is scaled to units of |cm^-1| and the parameters of
-the CylinderModel are the following:
-
-%(parameters)s
-
-The output of the 1D scattering intensity function for randomly oriented
-cylinders is then given by
-
-.. math::
-
-    P(q) = \frac{\text{scale}}{V}
-        \int_0^{\pi/2} f^2(q,\alpha) \sin \alpha d\alpha + \text{background}
-
-The *theta* and *phi* parameters are not used for the 1D output. Our
-implementation of the scattering kernel and the 1D scattering intensity
-use the c-library from NIST.
-
-Validation of the CylinderModel
--------------------------------
-
-Validation of our code was done by comparing the output of the 1D model
-to the output of the software provided by the NIST (Kline, 2006).
-Figure :num:`figure #CylinderModel-compare` shows a comparison of
-the 1D output of our model and the output of the NIST software.
-
-.. _CylinderModel-compare:
-
-.. figure:: img/image065.JPG
-
-    Comparison of the SasView scattering intensity for a cylinder with the
-    output of the NIST SANS analysis software.
-    The parameters were set to: *Scale* = 1.0, *Radius* = 20 |Ang|,
-    *Length* = 400 |Ang|, *Contrast* = 3e-6 |Ang^-2|, and
-    *Background* = 0.01 |cm^-1|.
-
-In general, averaging over a distribution of orientations is done by
-evaluating the following
-
-.. math::
-
-    P(q) = \int_0^{\pi/2} d\phi
-        \int_0^\pi p(\theta, \phi) P_0(q,\alpha) \sin \theta d\theta
-
-
-where $p(\theta,\phi)$ is the probability distribution for the orientation
-and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented
-system. Since we have no other software to compare the implementation of
-the intensity for fully oriented cylinders, we can compare the result of
-averaging our 2D output using a uniform distribution $p(\theta, \phi) = 1.0$.
-Figure :num:`figure #CylinderModel-crosscheck` shows the result of
-such a cross-check.
-
-.. _CylinderModel-crosscheck:
-
-.. figure:: img/image066.JPG
-
-    Comparison of the intensity for uniformly distributed cylinders
-    calculated from our 2D model and the intensity from the NIST SANS
-    analysis software.
-    The parameters used were: *Scale* = 1.0, *Radius* = 20 |Ang|,
-    *Length* = 400 |Ang|, *Contrast* = 3e-6 |Ang^-2|, and
-    *Background* = 0.0 |cm^-1|.
-"""
-
-from numpy import pi, inf
-
-name = "cylinder"
-title = "Cylinder with uniform scattering length density"
-description = """
-     f(q)= 2*(sldCyl - sldSolv)*V*sin(qLcos(alpha/2))
-            /[qLcos(alpha/2)]*J1(qRsin(alpha/2))/[qRsin(alpha)]
-
-            P(q,alpha)= scale/V*f(q)^(2)+background
-            V: Volume of the cylinder
-            R: Radius of the cylinder
-            L: Length of the cylinder
-            J1: The bessel function
-            alpha: angle betweenthe axis of the
-            cylinder and the q-vector for 1D
-            :the ouput is P(q)=scale/V*integral
-            from pi/2 to zero of...
-            f(q)^(2)*sin(alpha)*dalpha+ bkg
-    """
-
-parameters = [
-#   [ "name", "units", default, [lower, upper], "type",
-#     "description" ],
-    [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "",
-      "Cylinder scattering length density" ],
-    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "",
-      "Solvent scattering length density" ],
-    [ "radius", "Ang",  20, [0, inf], "volume",
-      "Cylinder radius" ],
-    [ "length", "Ang",  400, [0, inf], "volume",
-      "Cylinder length" ],
-    [ "theta", "degrees", 60, [-inf, inf], "orientation",
-      "In plane angle" ],
-    [ "phi", "degrees", 60, [-inf, inf], "orientation",
-      "Out of plane angle" ],
-    ]
-
-source = [ "lib/J1.c", "lib/gauss76.c", "lib/cylkernel.c" ]
-
-form_volume = """
-    return M_PI*radius*radius*length;
-    """
-
-Iq = """
-    const real h = REAL(0.5)*length;
-    real summ = REAL(0.0);
-    for (int i=0; i<76 ;i++) {
-        //const real zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
-        const real zi = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);
-        summ += Gauss76Wt[i] * CylKernel(q, radius, h, zi);
-    }
-    //const real form = (uplim-lolim)/2.0*summ;
-    const real form = summ * M_PI_4;
-
-    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
-    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
-    // The additional volume factor is for polydisperse volume normalization.
-    const real s = (sld - solvent_sld) * form_volume(radius, length);
-    return REAL(1.0e-4) * form * s * s;
-    """
-
-Iqxy = """
-    real sn, cn; // slots to hold sincos function output
-
-    // Compute angle alpha between q and the cylinder axis
-    SINCOS(theta*M_PI_180, sn, cn);
-    // # The following correction factor exists in sasview, but it can't be
-    // # right, so we are leaving it out for now.
-    // const real correction = fabs(cn)*M_PI_2;
-    const real q = sqrt(qx*qx+qy*qy);
-    const real cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
-    const real alpha = acos(cos_val);
-
-    // The following is CylKernel() / sin(alpha), but we are doing it in place
-    // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit
-    // faster since we don't mulitply and divide sin(alpha).
-    SINCOS(alpha, sn, cn);
-    const real besarg = q*radius*sn;
-    const real siarg = REAL(0.5)*q*length*cn;
-    // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
-    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg);
-    const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg);
-    const real form = REAL(4.0)*bj*bj*si*si;
-
-    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
-    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
-    // The additional volume factor is for polydisperse volume normalization.
-    const real s = (sld - solvent_sld) * form_volume(radius, length);
-    return REAL(1.0e-4) * form * s * s; // * correction;
-    """
-
-def ER(radius, length):
-    ddd = 0.75*radius*(2*radius*length + (length+radius)*(length+pi*radius))
-    return 0.5 * (ddd)**(1./3.)
-
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/sasview_model.py	(revision 87985ca5989a70a24a7e9ac09a971d544fe16dfe)
@@ -1,2 +1,21 @@
+"""
+Sasview model constructor.
+
+Given a module defining an OpenCL kernel such as sasmodels.models.cylinder,
+create a sasview model class to run that kernel as follows::
+
+    from sasmodels.sasview_model import make_class
+    from sasmodels.models import cylinder
+    CylinderModel = make_class(cylinder, dtype='single')
+
+The model parameters for sasmodels are different from those in sasview.
+When reloading previously saved models, the parameters should be converted
+using :func:`sasmodels.convert.convert`.
+"""
+
+# TODO: add a sasview=>sasmodels parameter translation layer
+# this will allow us to use the new sasmodels as drop in replacements, and
+# delay renaming parameters until all models have been converted.
+
 import math
 from copy import deepcopy
@@ -16,7 +35,6 @@
     """
     Load the sasview model defined in *kernel_module*.
-    :param kernel_module:
-    :param dtype:
-    :return:
+
+    Returns a class that can be used directly as a sasview model.
     """
     model =  load_model(kernel_module, dtype=dtype)
Index: setup_py2exe.py
===================================================================
--- setup_py2exe.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ setup_py2exe.py	(revision 87985ca5989a70a24a7e9ac09a971d544fe16dfe)
@@ -40,5 +40,5 @@
 import sys
 
-sys.dont_write_bytecode = True
+#sys.dont_write_bytecode = True
 
 # Force build before continuing
@@ -228,5 +228,5 @@
 # Specify required packages to bundle in the executable image.
 packages = ['numpy', 'scipy', 'matplotlib', 'pytz', 'pyparsing',
-            'periodictable', 'bumps.names', 'dream'
+            'periodictable', 'bumps', 'sasmodels', 'pyopencl',
             ]
 
