Index: sansmodels/src/sans/models/c_smearer/smearer2d_helper.cpp
===================================================================
--- sansmodels/src/sans/models/c_smearer/smearer2d_helper.cpp	(revision 87615a4840b7cfe8cb701bebbed7fd96effdc759)
+++ sansmodels/src/sans/models/c_smearer/smearer2d_helper.cpp	(revision 87615a4840b7cfe8cb701bebbed7fd96effdc759)
@@ -0,0 +1,77 @@
+/**
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	If you use DANSE applications to do scientific research that leads to
+	publication, we ask that you acknowledge the use of the software with the
+	following sentence:
+
+	"This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+	copyright 2009, University of Tennessee
+ */
+#include "smearer2d_helper.hh"
+#include <stdio.h>
+#include <math.h>
+using namespace std;
+
+/**
+ * Constructor for BaseSmearer
+ *
+ * binning
+ * @param qx: array of Qx values
+ * @param qy: array of Qy values
+ * @param nrbins: number of r bins
+ * @param nphibins: number of phi bins
+ */
+Smearer_helper :: Smearer_helper(int npoints, double* qx, double* qy,
+		double* dqx, double* dqy, double rlimit, int nrbins, int nphibins) {
+	// Number of bins
+	this->npoints = npoints;
+	this->rlimit = rlimit;
+	this->nrbins = nrbins;
+	this->nphibins = nphibins;
+	this->qx_values = qx;
+	this->qy_values = qy;
+	this->dqx_values = dqx;
+	this->dqy_values = dqy;
+};
+
+/**
+ * Compute the point smearing matrix
+ */
+void Smearer_helper :: smear2d(double *weights, double *qx_out, double *qy_out){
+
+	double rbin_size = rlimit / double(nrbins);
+	double phibin_size = 0.0;
+	int tot_nbins = nrbins * nphibins;
+	double rbin = 0.0;
+	double phibin = 0.0;
+	double qr = 0.0;
+	double qphi = 0.0;
+	double Pi = 4.0*atan(1.0);
+	// Loop over q-values and multiply apply matrix
+
+	for(int i=0; i<nrbins; i++){
+		rbin = rbin_size * (double(i) + 0.5);
+		for(int j=0; j<nphibins; j++){
+			phibin_size =  2.0 * Pi / double(nphibins);
+			phibin = phibin_size * (double(j));
+			for(int q_i=0; q_i<npoints; q_i++){
+				qr = sqrt(qx_values[q_i]*qx_values[q_i] + qy_values[q_i]*qy_values[q_i]);
+				qphi = atan(qy_values[q_i]/qx_values[q_i]);
+				qx_out[q_i + npoints*(nrbins * j + i)] = (rbin*dqx_values[q_i]*cos(phibin) + qr)*cos(qphi)-
+										rbin*dqy_values[q_i]*sin(phibin)*sin(qphi);
+				qy_out[q_i + npoints*(nrbins * j + i)] = (rbin*dqx_values[q_i]*cos(phibin) + qr)*sin(qphi)+
+										rbin*dqy_values[q_i]*sin(phibin)*cos(qphi);
+				if (q_i==0){
+					weights[nrbins * j + i] = exp(-0.5 * ((rbin - rbin_size / 2.0) *
+							(rbin - rbin_size / 2.0)))- exp(-0.5 * ((rbin + rbin_size / 2.0 ) *
+									(rbin + rbin_size / 2.0)));
+			}
+		}
+	}
+	}
+}
+
Index: sansmodels/src/sans/models/c_smearer/smearer2d_helper.hh
===================================================================
--- sansmodels/src/sans/models/c_smearer/smearer2d_helper.hh	(revision 87615a4840b7cfe8cb701bebbed7fd96effdc759)
+++ sansmodels/src/sans/models/c_smearer/smearer2d_helper.hh	(revision 87615a4840b7cfe8cb701bebbed7fd96effdc759)
@@ -0,0 +1,50 @@
+/**
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	If you use DANSE applications to do scientific research that leads to
+	publication, we ask that you acknowledge the use of the software with the
+	following sentence:
+
+	"This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+	copyright 2009, University of Tennessee
+ */
+#ifndef SMEARER2D_HELPER_CLASS_H
+#define SMEARER2D_HELPER_CLASS_H
+
+#include <vector>
+
+using namespace std;
+
+/**
+ * Base smearer class, implementing the matrix multiplication only
+ */
+class Smearer_helper {
+protected:
+	// Smearing matrix
+	vector<double>* weights;
+	vector<double>* qx_out;
+	vector<double>* qy_out;
+	// Q vector
+	double* qx_values;
+	double* qy_values;
+	double* dqx_values;
+	double* dqy_values;
+    // r limit
+    double rlimit;
+    // Numbers
+    int npoints;
+    int nrbins;
+    int nphibins;
+
+public:
+    // Constructor
+	Smearer_helper(int npoints, double* qx, double* qy, double* dqx, double* dqy,
+			double rlimit, int nrbins, int nphibins);
+	// Smear function
+	virtual void smear2d(double *weights, double *qx_out, double *qy_out);
+};
+
+#endif
Index: sansmodels/src/sans/models/c_smearer/smearer2d_helper_module.cpp
===================================================================
--- sansmodels/src/sans/models/c_smearer/smearer2d_helper_module.cpp	(revision 87615a4840b7cfe8cb701bebbed7fd96effdc759)
+++ sansmodels/src/sans/models/c_smearer/smearer2d_helper_module.cpp	(revision 87615a4840b7cfe8cb701bebbed7fd96effdc759)
@@ -0,0 +1,130 @@
+/**
+ *  Smearer module to perform point and slit smearing calculations
+
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	See the license text in license.txt
+
+	copyright 2009, University of Tennessee
+ */
+#include <Python.h>
+#include <stdio.h>
+#include <smearer2d_helper.hh>
+
+// Utilities
+#define INVECTOR(obj,buf,len)										\
+    do { \
+        int err = PyObject_AsReadBuffer(obj, (const void **)(&buf), &len); \
+        if (err < 0) return NULL; \
+        len /= sizeof(*buf); \
+    } while (0)
+
+#define OUTVECTOR(obj,buf,len) \
+    do { \
+        int err = PyObject_AsWriteBuffer(obj, (void **)(&buf), &len); \
+        if (err < 0) return NULL; \
+        len /= sizeof(*buf); \
+    } while (0)
+
+
+/**
+ * Delete a Smearer2d object
+ */
+void del_smearer2d_helper(void *ptr){
+	Smearer_helper* smearer2d_helper = static_cast<Smearer_helper *>(ptr);
+	delete smearer2d_helper;
+	return;
+}
+
+/**
+/**
+ * Create a QSmearer as a python object by supplying a q array
+ */
+PyObject * new_Smearer_helper(PyObject *, PyObject *args) {
+	PyObject *qx_values_obj;
+	PyObject *qy_values_obj;
+	PyObject *dqx_values_obj;
+	PyObject *dqy_values_obj;
+	//Py_ssize_t n_q;
+	//PyObject rlimit_obj;
+	//PyObject npoints_obj;
+	//PyObject nrbins_obj;
+	//PyObject nphibins_obj;
+	double* qx_values;
+	double* qy_values;
+	double* dqx_values;
+	double* dqy_values;
+	double rlimit;
+	int npoints;
+	int nrbins;
+	int nphibins;
+
+
+	if (!PyArg_ParseTuple(args, "OOOOdiii", &qx_values_obj, &qy_values_obj, &dqx_values_obj, &dqy_values_obj, &rlimit, &nrbins, &nphibins, &npoints)) return NULL;
+	OUTVECTOR(qx_values_obj, qx_values, npoints);
+	OUTVECTOR(qy_values_obj, qy_values, npoints);
+	OUTVECTOR(dqx_values_obj, dqx_values, npoints);
+	OUTVECTOR(dqy_values_obj, dqy_values, npoints);
+	Smearer_helper* smearer2d_helper = new Smearer_helper(npoints,qx_values,qy_values,dqx_values,dqy_values, rlimit, nrbins, nphibins);
+	return PyCObject_FromVoidPtr(smearer2d_helper, del_smearer2d_helper);
+}
+
+/**
+ * Smear the given input according to a given smearer object
+ */
+PyObject * smear2d_input(PyObject *, PyObject *args) {
+	PyObject *weights_obj;
+	PyObject *qx_out_obj;
+	PyObject *qy_out_obj;
+	PyObject *smear_obj;
+	Py_ssize_t w_out;
+	Py_ssize_t n_out;
+	double *weights;
+	double *qx_out;
+	double *qy_out;
+
+
+
+	if (!PyArg_ParseTuple(args, "OOOO",  &smear_obj, &weights_obj, &qx_out_obj, &qy_out_obj)) return NULL;
+	OUTVECTOR(weights_obj, weights, w_out);
+	OUTVECTOR(qx_out_obj, qx_out, n_out);
+	OUTVECTOR(qy_out_obj, qy_out, n_out);
+
+	// Sanity check
+	//if(n_in!=n_out) return Py_BuildValue("i",-1);
+
+	// Set the array pointers
+	void *temp = PyCObject_AsVoidPtr(smear_obj);
+	Smearer_helper* s = static_cast<Smearer_helper *>(temp);
+
+	s->smear2d(weights, qx_out, qy_out);
+	//return PyCObject_FromVoidPtr(s, del_smearer2d_helper);
+	return Py_BuildValue("i",1);
+}
+
+
+/**
+ * Define module methods
+ */
+static PyMethodDef module_methods[] = {
+	{"new_Smearer_helper", (PyCFunction)new_Smearer_helper, METH_VARARGS,
+		  "Create a new Q Smearer object"},
+	{"smearer2d_helper",(PyCFunction)smear2d_input, METH_VARARGS,
+		  "Smear the given input array"},
+    {NULL}
+};
+
+
+#ifndef PyMODINIT_FUNC	/* declarations for DLL import/export */
+#define PyMODINIT_FUNC void
+#endif
+PyMODINIT_FUNC
+initsmearer2d_helper(void)
+{
+    PyObject* m;
+
+    m = Py_InitModule3("smearer2d_helper", module_methods,
+                       "Smearer2d_helper module");
+}
Index: sansmodels/src/sans/models/smearing_2d.py
===================================================================
--- sansmodels/src/sans/models/smearing_2d.py	(revision 6e8b436b9448bd5d8b9b0fc4d82be7766f274b7e)
+++ sansmodels/src/sans/models/smearing_2d.py	(revision 87615a4840b7cfe8cb701bebbed7fd96effdc759)
@@ -10,4 +10,5 @@
 import numpy
 import math
+
 ## Singular point
 SIGMA_ZERO = 1.0e-010
@@ -16,6 +17,6 @@
 LIMIT = 3.0
 ## Defaults
-R_BIN = {'Xhigh':10.0, 'High':5.0,'Med':5.0,'Low':3.0}
-PHI_BIN ={'Xhigh':20.0,'High':12.0,'Med':6.0,'Low':4.0}   
+R_BIN = {'Xhigh':10, 'High':5,'Med':5,'Low':3}
+PHI_BIN ={'Xhigh':20,'High':12,'Med':6,'Low':4}   
 
 class Smearer2D:
@@ -25,5 +26,5 @@
      
     def __init__(self, data=None, model=None, index=None, 
-                 limit=LIMIT, accuracy='Low', coords='polar'):
+                 limit=LIMIT, accuracy='Low', coords='polar', engine='python'):
         """
         Assumption: equally spaced bins in dq_r, dq_phi space.
@@ -35,4 +36,5 @@
         :param nphi: number of bins in dq_phi-axis
         :param coord: coordinates [string], 'polar' or 'cartesian'
+        :param engine: engine name [string]; 'c' or 'numpy'
         """
         ## data
@@ -51,4 +53,5 @@
         self.coords = coords
         self.smearer = True
+        self._engine = engine
         
         
@@ -119,35 +122,5 @@
         Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, 
         then find smeared intensity
-        For the default values, this is equivalent (but by using numpy array 
-        the speed optimized by a factor of ten)to the following: ::
-       
-       
-        Remove the singular points if exists
-        self.dqx_data[self.dqx_data==0]=SIGMA_ZERO
-        self.dqy_data[self.dqy_data==0]=SIGMA_ZERO
-        
-        for phi in range(0,4):
-            for r in range(0,5):
-                n = (phi)*5+(r)
-                r = r+0.25
-                dphi = phi*2.0*math.pi/4.0 + numpy.arctan( \
-                        self.qy_data[index_model]/self.dqy_data[index_model]/ \
-                        self.qx_data[index_model]*/self.dqx_data[index_model])
-                dq = r*sqrt( self.dqx_data[index_model]*\
-                        self.dqx_data[index_model] \
-                    + self.dqy_data[index_model]*self.dqy_data[index_model] )
-                #integrant of exp(-0.5*r*r) r dr at each bins : 
-                The integration may not need.
-                weight_res[n] = e^{(-0.5*((r-0.25)*(r-0.25)))}- \
-                                e^{(-0.5*((r-0.25)*(r-0.25)))}
-                #if phi != 0 and r != 0:
-                qx_res = numpy.append(qx_res,self.qx_data[index_model]+ \
-                                                            dq * cos(dphi))
-                qy_res = numpy.append(qy_res,self.qy_data[index_model]+ \
-                                                            dq * sin(dphi))
-                    
-        Then compute I(qx_res,qy_res) and do weighted averaging. 
-        
-        """
+        """    
         valid = self.get_data()
         if valid == None:
@@ -163,9 +136,4 @@
         nr = self.nr[self.accuracy]
         nphi = self.nphi[self.accuracy]
-
-        # data length in the range of self.index
-        len_data = len(self.qx_data)
-        len_datay = len(self.qy_data)
-
         # Number of bins in the dqr direction (polar coordinate of dqx and dqy)
         bin_size = self.limit / nr
@@ -173,68 +141,82 @@
         # in dq_r-direction times # of bins in dq_phi-direction
         n_bins = nr * nphi
-        # Mean values of dqr at each bins ,starting from the half of bin size
-        r = bin_size / 2.0 + numpy.arange(nr) * bin_size
-        # mean values of qphi at each bines
-        phi = numpy.arange(nphi)
-        dphi = phi * 2.0 * math.pi / nphi
-        dphi = dphi.repeat(nr)
-
-        ## Transform to polar coordinate, 
-        #  and set dphi at each data points ; 1d array
-        dphi = dphi.repeat(len_data)
-        q_phi = self.qy_data / self.qx_data
-        
-        # Starting angle is different between polar and cartesian coordinates.
-        #if self.coords != 'polar':
-        #    dphi += numpy.arctan( q_phi * self.dqx_data/ \
-        #                     self.dqy_data).repeat(n_bins).reshape(len_data,\
-        #                                        n_bins).transpose().flatten()
-
-        # The angle (phi) of the original q point
-        q_phi = numpy.arctan(q_phi).repeat(n_bins).reshape(len_data,\
-                                                n_bins).transpose().flatten()
-        ## Find Gaussian weight for each dq bins: The weight depends only 
-        #  on r-direction (The integration may not need)
-        weight_res = numpy.exp(-0.5 * ((r - bin_size / 2.0) * \
-                                    (r - bin_size / 2.0)))- \
-                                    numpy.exp(-0.5 * ((r + bin_size / 2.0 ) *\
-                                    (r + bin_size / 2.0)))
-        # No needs of normalization here.
-        #weight_res /= numpy.sum(weight_res)
-        weight_res = weight_res.repeat(nphi).reshape(nr, nphi)
-
-        weight_res = weight_res.transpose().flatten()
-        
-        ## Set dr for all dq bins for averaging
-        dr = r.repeat(nphi).reshape(nr,nphi).transpose().flatten()
-        ## Set dqr for all data points
-        dqx = numpy.outer(dr,self.dqx_data).flatten()
-        dqy = numpy.outer(dr,self.dqy_data).flatten()
-
-        qx = self.qx_data.repeat(n_bins).reshape(len_data,\
-                                             n_bins).transpose().flatten()
-        qy = self.qy_data.repeat(n_bins).reshape(len_data,\
-                                             n_bins).transpose().flatten()
-
-        # The polar needs rotation by -q_phi
-        if self.coords == 'polar':
-            q_r = numpy.sqrt(qx * qx + qy * qy)
-            qx_res = ((dqx*numpy.cos(dphi) + q_r) * numpy.cos(-q_phi) +\
-                           dqy*numpy.sin(dphi) * numpy.sin(-q_phi))
-            qy_res = (-(dqx*numpy.cos(dphi) + q_r) * numpy.sin(-q_phi) +\
-                           dqy*numpy.sin(dphi) * numpy.cos(-q_phi))
+        # data length in the range of self.index
+        len_data = len(self.qx_data)
+        len_datay = len(self.qy_data)
+
+        if self._engine == 'c' and self.coords == 'polar':
+            try:
+                import sans_extension.smearer2d_helper as smearer2dc
+                smearc = smearer2dc.new_Smearer_helper(self.qx_data, self.qy_data,
+                                              self.dqx_data, self.dqy_data,
+                                              self.limit, nr, nphi, int(len_data))
+                weight_res = numpy.zeros(nr * nphi )
+                qx_res = numpy.zeros(nr * nphi * int(len_data))
+                qy_res = numpy.zeros(nr * nphi * int(len_data))
+                smearer2dc.smearer2d_helper(smearc,weight_res, qx_res, qy_res)
+            except:
+                raise 
         else:
-            qx_res = qx +  dqx*numpy.cos(dphi)
-            qy_res = qy +  dqy*numpy.sin(dphi)
-
+            # Mean values of dqr at each bins ,starting from the half of bin size
+            r = bin_size / 2.0 + numpy.arange(nr) * bin_size
+            # mean values of qphi at each bines
+            phi = numpy.arange(nphi)
+            dphi = phi * 2.0 * math.pi / nphi
+            dphi = dphi.repeat(nr)
+    
+            ## Transform to polar coordinate, 
+            #  and set dphi at each data points ; 1d array
+            dphi = dphi.repeat(len_data)
+            q_phi = self.qy_data / self.qx_data
+            
+            # Starting angle is different between polar and cartesian coordinates.
+            #if self.coords != 'polar':
+            #    dphi += numpy.arctan( q_phi * self.dqx_data/ \
+            #                     self.dqy_data).repeat(n_bins).reshape(len_data,\
+            #                                        n_bins).transpose().flatten()
+    
+            # The angle (phi) of the original q point
+            q_phi = numpy.arctan(q_phi).repeat(n_bins).reshape(len_data,\
+                                                    n_bins).transpose().flatten()
+            ## Find Gaussian weight for each dq bins: The weight depends only 
+            #  on r-direction (The integration may not need)
+            weight_res = numpy.exp(-0.5 * ((r - bin_size / 2.0) * \
+                                        (r - bin_size / 2.0)))- \
+                                        numpy.exp(-0.5 * ((r + bin_size / 2.0 ) *\
+                                        (r + bin_size / 2.0)))
+            # No needs of normalization here.
+            #weight_res /= numpy.sum(weight_res)
+            weight_res = weight_res.repeat(nphi).reshape(nr, nphi)
+    
+            weight_res = weight_res.transpose().flatten()
+            
+            ## Set dr for all dq bins for averaging
+            dr = r.repeat(nphi).reshape(nr,nphi).transpose().flatten()
+            ## Set dqr for all data points
+            dqx = numpy.outer(dr,self.dqx_data).flatten()
+            dqy = numpy.outer(dr,self.dqy_data).flatten()
+    
+            qx = self.qx_data.repeat(n_bins).reshape(len_data,\
+                                                 n_bins).transpose().flatten()
+            qy = self.qy_data.repeat(n_bins).reshape(len_data,\
+                                                 n_bins).transpose().flatten()
+    
+            # The polar needs rotation by -q_phi
+            if self.coords == 'polar':
+                q_r = numpy.sqrt(qx * qx + qy * qy)
+                qx_res = ((dqx*numpy.cos(dphi) + q_r) * numpy.cos(-q_phi) +\
+                               dqy*numpy.sin(dphi) * numpy.sin(-q_phi))
+                qy_res = (-(dqx*numpy.cos(dphi) + q_r) * numpy.sin(-q_phi) +\
+                               dqy*numpy.sin(dphi) * numpy.cos(-q_phi))
+            else:
+                qx_res = qx +  dqx*numpy.cos(dphi)
+                qy_res = qy +  dqy*numpy.sin(dphi)
+            
         ## Evaluate all points
         val = self.model.evalDistribution([qx_res, qy_res]) 
-
         ## Reshape into 2d array to use numpy weighted averaging
         value_res= val.reshape(n_bins,len(self.qx_data))
-
         ## Averaging with Gaussian weighting: normalization included.
         value =numpy.average(value_res,axis=0, weights=weight_res)
-
         ## Return the smeared values in the range of self.index
         return value
