Index: sansmodels/src/sans/models/c_smearer/smearer.cpp
===================================================================
--- sansmodels/src/sans/models/c_smearer/smearer.cpp	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
+++ sansmodels/src/sans/models/c_smearer/smearer.cpp	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
@@ -0,0 +1,374 @@
+/**
+	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 "smearer.hh"
+#include <stdio.h>
+#include <math.h>
+using namespace std;
+
+
+/**
+ * Constructor for BaseSmearer
+ *
+ * @param qmin: minimum Q value
+ * @param qmax: maximum Q value
+ * @param nbins: number of Q bins
+ */
+BaseSmearer :: BaseSmearer(double qmin, double qmax, int nbins) {
+	// Number of bins
+	this->nbins = nbins;
+	this->qmin = qmin;
+	this->qmax = qmax;
+	// Flag to keep track of whether we have a smearing matrix or
+	// whether we need to compute one
+	has_matrix = false;
+	even_binning = true;
+};
+
+/**
+ * Constructor for BaseSmearer
+ *
+ * Used for uneven binning
+ * @param q: array of Q values
+ * @param nbins: number of Q bins
+ */
+BaseSmearer :: BaseSmearer(double* q, int nbins) {
+	// Number of bins
+	this->nbins = nbins;
+	this->q_values = q;
+	// Flag to keep track of whether we have a smearing matrix or
+	// whether we need to compute one
+	has_matrix = false;
+	even_binning = false;
+};
+
+/**
+ * Constructor for SlitSmearer
+ *
+ * @param width: slit width in Q units
+ * @param height: slit height in Q units
+ * @param qmin: minimum Q value
+ * @param qmax: maximum Q value
+ * @param nbins: number of Q bins
+ */
+SlitSmearer :: SlitSmearer(double width, double height, double qmin, double qmax, int nbins) :
+	BaseSmearer(qmin, qmax, nbins){
+	this->height = height;
+	this->width = width;
+};
+
+/**
+ * Constructor for SlitSmearer
+ *
+ * @param width: slit width in Q units
+ * @param height: slit height in Q units
+ * @param q: array of Q values
+ * @param nbins: number of Q bins
+ */
+SlitSmearer :: SlitSmearer(double width, double height, double* q, int nbins) :
+	BaseSmearer(q, nbins){
+	this->height = height;
+	this->width = width;
+};
+
+/**
+ * Constructor for QSmearer
+ *
+ * @param width: array slit widths for each Q point, in Q units
+ * @param qmin: minimum Q value
+ * @param qmax: maximum Q value
+ * @param nbins: number of Q bins
+ */
+QSmearer :: QSmearer(double* width, double qmin, double qmax, int nbins) :
+	BaseSmearer(qmin, qmax, nbins){
+	this->width = width;
+};
+
+/**
+ * Constructor for QSmearer
+ *
+ * @param width: array slit widths for each Q point, in Q units
+ * @param q: array of Q values
+ * @param nbins: number of Q bins
+ */
+QSmearer :: QSmearer(double* width, double* q, int nbins) :
+	BaseSmearer(q, nbins){
+	this->width = width;
+};
+
+/**
+ * Compute the slit smearing matrix
+ *
+ * For even binning (q_min to q_max with nbins):
+ *
+ *   step = (q_max-q_min)/(nbins-1)
+ *   first bin goes from q_min to q_min+step
+ *   last bin goes from q_max to q_max+step
+ *
+ * For binning according to q array:
+ *
+ * Each q point represents a bin going from half the distance between it
+ * and the previous point to half the distance between it and the next point.
+ *
+ *    Example: bin i goes from (q_values[i-1]+q_values[i])/2 to (q_values[i]+q_values[i+1])/2
+ *
+ * The exceptions are the first and last bins, which are centered at the first and
+ * last q-values, respectively. The width of the first and last bins is the distance between
+ * their respective neighboring q-value.
+ */
+void SlitSmearer :: compute_matrix(){
+
+	weights = new vector<double>(nbins*nbins,0);
+
+	// Check the length of the data
+	if (nbins<2) return;
+	int npts_h = height>0.0 ? npts : 1;
+	int npts_w = width>0.0 ? npts : 1;
+
+	// If both height and width are great than zero,
+	// modify the number of points in each direction so
+	// that the total number of points is still what
+	// the user would expect (downgrade resolution)
+	//if(npts_h>1 && npts_w>1){
+	//	npts_h = (int)ceil(sqrt((double)npts));
+	//	npts_w = npts_h;
+	//}
+	double shift_h, shift_w, hbin_size, wbin_size;
+	// Make sure height and width are all positive (FWMH/2)
+	// Assumption; height and width are all same for all q points
+	if(npts_h == 1){
+		shift_h = 0.0;
+	} else {
+		shift_h = fabs(height);
+	}
+	if(npts_w == 1){
+		shift_w = 0.0;
+	} else {
+		shift_w = fabs(width);
+	}
+	// size of the h bin and w bin
+	hbin_size = shift_h / nbins;
+	wbin_size = shift_w / nbins;
+
+	// Loop over all q-values
+	for(int i=0; i<nbins; i++) {
+		// Find Weights
+		// Find q where the resolution smearing calculation of I(q) occurs
+		double q, q_min, q_max, q_0;
+		get_bin_range(i, &q, &q_min, &q_max);
+		// Block q becomes <=0
+		if (q <= 0){
+			continue;
+		}
+		bool last_qpoint = true;
+		// Find q[0] value to normalize the weight later,
+		//  otherwise, we will have a precision problem.
+		if (i == 0){
+			q_0 = q;
+		}
+		// Loop over all qj-values
+		bool first_w = true;
+		for(int j=0; j<nbins; j++) {
+			double q_j, q_high, q_low;
+			// Calculate bin size of q_j
+			get_bin_range(j, &q_j, &q_low, &q_high);
+			// Block q_j becomes <=0
+			if (q_j <= 0){
+				continue;
+			}
+			// Check q_low that can not be negative.
+			if (q_low < 0.0){
+				q_low = 0.0;
+			}
+			// default parameter values
+			(*weights)[i*nbins+j] = 0.0;
+			// protect for negative q
+			if (q <= 0.0 || q_j <= 0.0){
+					continue;
+			}
+			double shift_w = 0.0;
+			// Condition: zero slit smear.
+			if (npts_w == 1 && npts_h == 1){
+				if(q_j == q) {
+					(*weights)[i*nbins+j] = 1.0;
+				}
+			}
+			//Condition:Smear weight integration for width >0 when the height (=0) does not present.
+			//Or height << width.
+			else if((npts_w!=1 && npts_h == 1)|| (npts_w!=1 && npts_h != 1 && width/height > 100.0)){
+				shift_w = width;
+				//del_w = width/((double)npts_w-1.0);
+				double q_shifted_low = q - shift_w;
+				// High limit of the resolution range
+				double q_shifted_high = q + shift_w;
+				// Go through all the q_js for weighting those points
+				if(q_j >= q_shifted_low && q_j <= q_shifted_high) {
+					// The weighting factor comes,
+					// Give some weight (delq_bin) for the q_j within the resolution range
+					// Weight should be same for all qs except
+					// for the q bin size at j.
+					// Note that the division by q_0 is only due to the precision problem
+					// where q_high - q_low gets to very small.
+					// Later, it will be normalized again.
+					(*weights)[i*nbins+j] += (q_high - q_low)/q_0 ;
+				}
+			}
+			else{
+				// Loop for width (;Height is analytical.)
+				// Condition: height >>> width, otherwise, below is not accurate enough.
+				// Smear weight numerical iteration for width >0 when the height (>0) presents.
+				// When width = 0, the numerical iteration will be skipped.
+				// The resolution calculation for the height is done by direct integration,
+				// assuming the I(q'=sqrt(q_j^2-(q+shift_w)^2)) is constant within a q' bin, [q_high, q_low].
+				// In general, this weight numerical iteration for width >0 might be a rough approximation,
+				// but it must be good enough when height >>> width.
+				for(int k=(-npts_w + 1); k<npts_w; k++){
+					if(npts_w!=1){
+						shift_w = width/((double)npts_w-1.0)*(double)k;
+					}
+					// For each q-value, compute the weight of each other q-bin
+					// in the I(q) array
+					// Low limit of the resolution range
+					double q_shift = q + shift_w;
+					if (q_shift < 0.0){
+						q_shift = 0.0;
+					}
+					double q_shifted_low = q_shift;
+					// High limit of the resolution range
+					double q_shifted_high = sqrt(q_shift * q_shift + shift_h * shift_h);
+
+
+					// Go through all the q_js for weighting those points
+					if(q_j >= q_shifted_low && q_j <= q_shifted_high) {
+						// The weighting factor comes,
+						// Give some weight (delq_bin) for the q_j within the resolution range
+						// Weight should be same for all qs except
+						// for the q bin size at j.
+						// Note that the division by q_0 is only due to the precision problem
+						// where q_high - q_low gets to very small.
+						// Later, it will be normalized again.
+
+						double q_shift_min = q - width;
+
+						double u = (q_j * q_j - (q_shift) * (q_shift));
+						// The fabs below are not necessary but in case: the weight should never be imaginary.
+						// At the edge of each sub_width. weight += u(at q_high bin) - u(0), where u(0) = 0,
+						// and weighted by (2.0* npts_w -1.0)once for each q.
+						if (q == q_j) {
+							if (k==0)
+								(*weights)[i*nbins+j] += (sqrt(fabs((q_high)*(q_high)-q_shift * q_shift)))/q_0 * (2.0*double(npts_w)-1.0);
+						}
+						// For the rest of sub_width. weight += u(at q_high bin) - u(at q_low bin)
+						else if (u > 0.0){
+							(*weights)[i*nbins+j] += (sqrt(fabs((q_high)*(q_high)- q_shift * q_shift))-sqrt(fabs((q_low)*(q_low)- q_shift * q_shift)))/q_0 ;
+						}
+					}
+				}
+			}
+		}
+	}
+};
+
+/**
+ * Compute the point smearing matrix
+ */
+void QSmearer :: compute_matrix(){
+	weights = new vector<double>(nbins*nbins,0);
+
+	// Loop over all q-values
+	double step = (qmax-qmin)/((double)nbins-1.0);
+	double q, q_min, q_max;
+	double q_j, q_jmax, q_jmin;
+	for(int i=0; i<nbins; i++) {
+		get_bin_range(i, &q, &q_min, &q_max);
+
+		for(int j=0; j<nbins; j++) {
+			get_bin_range(j, &q_j, &q_jmin, &q_jmax);
+
+			// Compute the fraction of the Gaussian contributing
+			// to the q_j bin between q_jmin and q_jmax
+			double value =  erf( (q_jmax-q)/(sqrt(2.0)*width[i]) );
+        	value -= erf( (q_jmin-q)/(sqrt(2.0)*width[i]) );
+        	(*weights)[i*nbins+j] += value;
+		}
+	}
+}
+
+/**
+ * Computes the Q range of a given bin of the Q distribution.
+ * The range is computed according the the data distribution that
+ * was given to the object at initialization.
+ *
+ * @param i: number of the bin in the distribution
+ * @param q: q-value of bin i
+ * @param q_min: lower bound of the bin
+ * @param q_max: higher bound of the bin
+ *
+ */
+int BaseSmearer :: get_bin_range(int i, double* q, double* q_min, double* q_max) {
+	if (even_binning) {
+		double step = (qmax-qmin)/((double)nbins-1.0);
+		*q = qmin + (double)i*step;
+		*q_min = *q - 0.5*step;
+		*q_max = *q + 0.5*step;
+		return 1;
+	} else if (i>=0 && i<nbins) {
+		*q = q_values[i];
+		if (i==0) {
+			double step = (q_values[1]-q_values[0])/2.0;
+			*q_min = *q - step;
+			*q_max = *q + step;
+		} else if (i==nbins-1) {
+			double step = (q_values[i]-q_values[i-1])/2.0;
+			*q_min = *q - step;
+			*q_max = *q + step;
+		} else {
+			*q_min = *q - (q_values[i]-q_values[i-1])/2.0;
+			*q_max = *q + (q_values[i+1]-q_values[i])/2.0;
+		}
+		return 1;
+	}
+	return -1;
+}
+
+/**
+ * Perform smearing by applying the smearing matrix to the input Q array
+ */
+void BaseSmearer :: smear(double *iq_in, double *iq_out, int first_bin, int last_bin){
+
+	// If we haven't computed the smearing matrix, do it now
+	if(!has_matrix) {
+		compute_matrix();
+		has_matrix = true;
+	}
+
+	// Loop over q-values and multiply apply matrix
+	for(int q_i=first_bin; q_i<=last_bin; q_i++){
+		double sum = 0.0;
+		double counts = 0.0;
+
+		for(int i=first_bin; i<=last_bin; i++){
+			// Skip if weight is less than 1e-03(this value is much smaller than
+			// the weight at the 3*sigma distance
+			// Will speed up a little bit...
+			if ((*weights)[q_i*nbins+i] < 1.0e-003){
+				continue;
+			}
+			sum += iq_in[i] * (*weights)[q_i*nbins+i];
+			counts += (*weights)[q_i*nbins+i];
+		}
+
+		// Normalize counts
+		iq_out[q_i] = (counts>0.0) ? sum/counts : 0.0;
+	}
+}
Index: sansmodels/src/sans/models/c_smearer/smearer.hh
===================================================================
--- sansmodels/src/sans/models/c_smearer/smearer.hh	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
+++ sansmodels/src/sans/models/c_smearer/smearer.hh	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
@@ -0,0 +1,97 @@
+/**
+	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 SMEARER_CLASS_H
+#define SMEARER_CLASS_H
+
+#include <vector>
+
+using namespace std;
+
+/**
+ * Base smearer class, implementing the matrix multiplication only
+ */
+class BaseSmearer {
+protected:
+	// Internal flag: true when the weights vector is filled
+	bool has_matrix;
+	// True when only qmin, qmax and nbins are given. Bins are equidistant
+	bool even_binning;
+	// Smearing matrix
+	vector<double>* weights;
+	// Q vector
+	double* q_values;
+    // Q_min (Min Q-value for I(q))
+    double qmin;
+    // Q_max (Max Q_value for I(q))
+    double qmax;
+    // Number of Q bins
+    int nbins;
+
+public:
+    // Constructor
+    BaseSmearer(double qmin, double qmax, int nbins);
+    BaseSmearer(double* q, int nbins);
+	// Smear function
+	virtual void smear(double *, double *, int, int);
+	// Compute the smearing matrix
+	virtual void compute_matrix(){};
+	// Utility function to check the number of bins
+	int get_nbins() { return nbins; }
+	// Get the q range of a particular bin
+	virtual int get_bin_range(int, double*, double*, double*);
+};
+
+
+/**
+ * Slit smearer class
+ */
+class SlitSmearer : public BaseSmearer {
+
+protected:
+    // Number of points used in the smearing computation
+    static const int npts   = 500;
+
+public:
+    // Slit width in Q units
+    double width;
+    // Slit height in Q units
+    double height;
+
+    // Constructor
+	SlitSmearer(double width, double height, double qmin, double qmax, int nbins);
+	SlitSmearer(double width, double height, double* q, int nbins);
+	// Compute the smearing matrix
+	virtual void compute_matrix();
+};
+
+
+/**
+ * Point smearer class
+ */
+class QSmearer : public BaseSmearer {
+
+protected:
+    // Standard deviation in Q [A-1]
+    double* width;
+
+public:
+
+    // Constructor
+    QSmearer(double* width, double qmin, double qmax, int nbins);
+    QSmearer(double* width, double* q, int nbins);
+	// Compute the smearing matrix
+	virtual void compute_matrix();
+};
+
+#endif
Index: sansmodels/src/sans/models/c_smearer/smearer_module.cpp
===================================================================
--- sansmodels/src/sans/models/c_smearer/smearer_module.cpp	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
+++ sansmodels/src/sans/models/c_smearer/smearer_module.cpp	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
@@ -0,0 +1,205 @@
+/**
+ *  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 <smearer.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 SlitSmearer object
+ */
+void del_slit_smearer(void *ptr){
+	SlitSmearer* smearer = static_cast<SlitSmearer *>(ptr);
+	delete smearer;
+	return;
+}
+
+/**
+ * Create a SlitSmearer as a python object
+ */
+PyObject * new_slit_smearer(PyObject *, PyObject *args) {
+	double width;
+	double height;
+	double qmin;
+	double qmax;
+	int nbins;
+
+	if (!PyArg_ParseTuple(args, "ddddi", &width, &height, &qmin, &qmax, &nbins)) return NULL;
+	SlitSmearer* smearer = new SlitSmearer(width, height, qmin, qmax, nbins);
+	return PyCObject_FromVoidPtr(smearer, del_slit_smearer);
+}
+
+/**
+ * Create a SlitSmearer as a python object by supplying a q array
+ */
+PyObject * new_slit_smearer_with_q(PyObject *, PyObject *args) {
+	PyObject *qvalues_obj;
+	Py_ssize_t n_q;
+	double* qvalues;
+	double width;
+	double height;
+
+	if (!PyArg_ParseTuple(args, "ddO", &width, &height, &qvalues_obj)) return NULL;
+	OUTVECTOR(qvalues_obj, qvalues, n_q);
+	SlitSmearer* smearer = new SlitSmearer(width, height, qvalues, n_q);
+	return PyCObject_FromVoidPtr(smearer, del_slit_smearer);
+}
+
+/**
+ * Delete a QSmearer object
+ */
+void del_q_smearer(void *ptr){
+	QSmearer* smearer = static_cast<QSmearer *>(ptr);
+	delete smearer;
+	return;
+}
+
+/**
+ * Create a QSmearer as a python object
+ */
+PyObject * new_q_smearer(PyObject *, PyObject *args) {
+	PyObject *width_obj;
+	Py_ssize_t n_w;
+	double* width;
+	double qmin;
+	double qmax;
+	int nbins;
+
+	if (!PyArg_ParseTuple(args, "Oddi", &width_obj, &qmin, &qmax, &nbins)) return NULL;
+	OUTVECTOR(width_obj, width, n_w);
+	QSmearer* smearer = new QSmearer(width, qmin, qmax, nbins);
+	return PyCObject_FromVoidPtr(smearer, del_q_smearer);
+}
+
+/**
+ * Create a QSmearer as a python object by supplying a q array
+ */
+PyObject * new_q_smearer_with_q(PyObject *, PyObject *args) {
+	PyObject *width_obj;
+	Py_ssize_t n_w;
+	double* width;
+	PyObject *qvalues_obj;
+	Py_ssize_t n_q;
+	double* qvalues;
+
+	if (!PyArg_ParseTuple(args, "OO", &width_obj, &qvalues_obj)) return NULL;
+	OUTVECTOR(width_obj, width, n_w);
+	OUTVECTOR(qvalues_obj, qvalues, n_q);
+	QSmearer* smearer = new QSmearer(width, qvalues, n_q);
+	return PyCObject_FromVoidPtr(smearer, del_q_smearer);
+}
+
+/**
+ * Smear the given input according to a given smearer object
+ */
+PyObject * smear_input(PyObject *, PyObject *args) {
+	PyObject *input_iq_obj;
+	PyObject *output_iq_obj;
+	PyObject *smear_obj;
+	Py_ssize_t n_in;
+	Py_ssize_t n_out;
+	double *input_iq;
+	double *output_iq;
+	int first_bin;
+	int last_bin;
+
+
+	if (!PyArg_ParseTuple(args, "OOOii", &smear_obj, &input_iq_obj, &output_iq_obj, &first_bin, &last_bin)) return NULL;
+	OUTVECTOR(input_iq_obj, input_iq, n_in);
+	OUTVECTOR(output_iq_obj, output_iq, n_out);
+
+	// Sanity check
+	if(n_in!=n_out) return Py_BuildValue("i",-1);
+
+	// Set the array pointers
+	void *temp = PyCObject_AsVoidPtr(smear_obj);
+	BaseSmearer* s = static_cast<BaseSmearer *>(temp);
+
+	if(n_in!=s->get_nbins()) {
+		return Py_BuildValue("i",-2);
+	}
+
+	s->smear(input_iq, output_iq, first_bin, last_bin);
+
+	return Py_BuildValue("i",1);
+}
+
+/**
+ * Get the q-value of a given bin
+ */
+PyObject * get_q(PyObject *, PyObject *args) {
+	PyObject *smear_obj;
+	int bin;
+
+	if (!PyArg_ParseTuple(args, "Oi", &smear_obj, &bin)) return NULL;
+
+	// Set the array pointers
+	void *temp = PyCObject_AsVoidPtr(smear_obj);
+	BaseSmearer* s = static_cast<BaseSmearer *>(temp);
+
+	if(s->get_nbins()<=0 || s->get_nbins()<=bin) {
+		return NULL;
+	}
+
+	double q, q_min, q_max;
+	if (s->get_bin_range(bin, &q, &q_min, &q_max)<0) {
+		return NULL;
+	}
+	return Py_BuildValue("d", q);
+}
+
+
+/**
+ * Define module methods
+ */
+static PyMethodDef module_methods[] = {
+	{"new_slit_smearer", (PyCFunction)new_slit_smearer, METH_VARARGS,
+		  "Create a new Slit Smearer object"},
+	{"new_slit_smearer_with_q", (PyCFunction)new_slit_smearer_with_q, METH_VARARGS,
+		  "Create a new Slit Smearer object by supplying an array of Q values"},
+	{"new_q_smearer", (PyCFunction)new_q_smearer, METH_VARARGS,
+		  "Create a new Q Smearer object"},
+	{"new_q_smearer_with_q", (PyCFunction)new_q_smearer_with_q, METH_VARARGS,
+		  "Create a new Q Smearer object by supplying an array of Q values"},
+	{"smear",(PyCFunction)smear_input, METH_VARARGS,
+		  "Smear the given input array"},
+	{"get_q",(PyCFunction)get_q, METH_VARARGS,
+		  "Get the q-value of a given bin"},
+    {NULL}
+};
+
+
+#ifndef PyMODINIT_FUNC	/* declarations for DLL import/export */
+#define PyMODINIT_FUNC void
+#endif
+PyMODINIT_FUNC
+initsmearer(void)
+{
+    PyObject* m;
+
+    m = Py_InitModule3("smearer", module_methods,
+                       "Q Smearing module");
+}
Index: sansmodels/src/sans/models/qsmearing.py
===================================================================
--- sansmodels/src/sans/models/qsmearing.py	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
+++ sansmodels/src/sans/models/qsmearing.py	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
@@ -0,0 +1,592 @@
+
+#####################################################################
+#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 2008, University of Tennessee
+######################################################################
+import numpy
+import math
+import logging
+import sys
+import sans_extension.smearer as smearer 
+from sans.models.smearing_2d import Smearer2D
+
+def smear_selection(data1D, model = None):
+    """
+    Creates the right type of smearer according 
+    to the data.
+
+    The canSAS format has a rule that either
+    slit smearing data OR resolution smearing data
+    is available. 
+    
+    For the present purpose, we choose the one that
+    has none-zero data. If both slit and resolution
+    smearing arrays are filled with good data 
+    (which should not happen), then we choose the
+    resolution smearing data. 
+    
+    :param data1D: Data1D object
+    :param model: sans.model instance
+    """
+    # Sanity check. If we are not dealing with a SANS Data1D
+    # object, just return None
+    if  data1D.__class__.__name__ not in ['Data1D', 'Theory1D']:
+        if data1D == None:
+            return None
+        elif data1D.dqx_data == None or data1D.dqy_data == None:
+            return None
+        return Smearer2D(data1D)
+    
+    if  not hasattr(data1D, "dx") and not hasattr(data1D, "dxl")\
+         and not hasattr(data1D, "dxw"):
+        return None
+    
+    # Look for resolution smearing data
+    _found_resolution = False
+    if data1D.dx is not None and len(data1D.dx) == len(data1D.x):
+        
+        # Check that we have non-zero data
+        if data1D.dx[0] > 0.0:
+            _found_resolution = True
+            #print "_found_resolution",_found_resolution
+            #print "data1D.dx[0]",data1D.dx[0],data1D.dxl[0]
+    # If we found resolution smearing data, return a QSmearer
+    if _found_resolution == True:
+        return QSmearer(data1D, model)
+
+    # Look for slit smearing data
+    _found_slit = False
+    if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x) \
+        and data1D.dxw is not None and len(data1D.dxw) == len(data1D.x):
+        
+        # Check that we have non-zero data
+        if data1D.dxl[0] > 0.0 or data1D.dxw[0] > 0.0:
+            _found_slit = True
+        
+        # Sanity check: all data should be the same as a function of Q
+        for item in data1D.dxl:
+            if data1D.dxl[0] != item:
+                _found_resolution = False
+                break
+            
+        for item in data1D.dxw:
+            if data1D.dxw[0] != item:
+                _found_resolution = False
+                break
+    # If we found slit smearing data, return a slit smearer
+    if _found_slit == True:
+        return SlitSmearer(data1D, model)
+    return None
+            
+
+class _BaseSmearer(object):
+    
+    def __init__(self):
+        self.nbins = 0
+        self.nbins_low = 0
+        self.nbins_high = 0
+        self._weights = None
+        ## Internal flag to keep track of C++ smearer initialization
+        self._init_complete = False
+        self._smearer = None
+        self.model = None
+        
+    def __deepcopy__(self, memo={}):
+        """
+        Return a valid copy of self.
+        Avoid copying the _smearer C object and force a matrix recompute
+        when the copy is used.  
+        """
+        result = _BaseSmearer()
+        result.nbins = self.nbins
+        return result
+
+    def _compute_matrix(self):
+        """
+        """
+        return NotImplemented
+
+    def get_bin_range(self, q_min=None, q_max=None):
+        """
+        
+        :param q_min: minimum q-value to smear
+        :param q_max: maximum q-value to smear
+         
+        """
+        # If this is the first time we call for smearing,
+        # initialize the C++ smearer object first
+        if not self._init_complete:
+            self._initialize_smearer()
+        if q_min == None:
+            q_min = self.min
+        if q_max == None:
+            q_max = self.max
+
+        _qmin_unsmeared, _qmax_unsmeared = self.get_unsmeared_range(q_min,
+                                                                     q_max)
+        _first_bin = None
+        _last_bin  = None
+
+        #step = (self.max - self.min) / (self.nbins - 1.0)
+        # Find the first and last bin number in all extrapolated and real data
+        try:
+            for i in range(self.nbins):
+                q_i = smearer.get_q(self._smearer, i)
+                if (q_i >= _qmin_unsmeared) and (q_i <= _qmax_unsmeared):
+                    # Identify first and last bin
+                    if _first_bin is None:
+                        _first_bin = i
+                    else:
+                        _last_bin  = i
+        except:
+            msg = "_BaseSmearer.get_bin_range: "
+            msg += " error getting range\n  %s" % sys.exc_value
+            raise RuntimeError, msg
+   
+        #  Find the first and last bin number only in the real data
+        _first_bin, _last_bin = self._get_unextrapolated_bin( \
+                                        _first_bin, _last_bin)
+
+        return _first_bin, _last_bin
+
+    def __call__(self, iq_in, first_bin = 0, last_bin = None):
+        """
+        Perform smearing
+        """
+        # If this is the first time we call for smearing,
+        # initialize the C++ smearer object first
+        if not self._init_complete:
+            self._initialize_smearer()
+
+        if last_bin is None or last_bin >= len(iq_in):
+            last_bin = len(iq_in) - 1
+        # Check that the first bin is positive
+        if first_bin < 0:
+            first_bin = 0
+
+        # With a model given, compute I for the extrapolated points and append
+        # to the iq_in
+        iq_in_temp = iq_in
+        if self.model != None:
+            temp_first, temp_last = self._get_extrapolated_bin( \
+                                                        first_bin, last_bin)
+            if self.nbins_low > 0:
+                iq_in_low = self.model.evalDistribution( \
+                                    numpy.fabs(self.qvalues[0:self.nbins_low]))
+            iq_in_high = self.model.evalDistribution( \
+                                            self.qvalues[(len(self.qvalues) - \
+                                            self.nbins_high - 1):])
+            # Todo: find out who is sending iq[last_poin] = 0.
+            if iq_in[len(iq_in) - 1] == 0:
+                iq_in[len(iq_in) - 1] = iq_in_high[0]
+            # Append the extrapolated points to the data points
+            if self.nbins_low > 0:                             
+                iq_in_temp = numpy.append(iq_in_low, iq_in)
+            if self.nbins_high > 0:
+                iq_in_temp = numpy.append(iq_in_temp, iq_in_high[1:])
+        else:
+            temp_first = first_bin
+            temp_last = last_bin
+            #iq_in_temp = iq_in
+
+        # Sanity check
+        if len(iq_in_temp) != self.nbins:
+            msg = "Invalid I(q) vector: inconsistent array "
+            msg += " length %d != %s" % (len(iq_in_temp), str(self.nbins))
+            raise RuntimeError, msg
+
+        # Storage for smeared I(q)   
+        iq_out = numpy.zeros(self.nbins)
+
+        smear_output = smearer.smear(self._smearer, iq_in_temp, iq_out,
+                                      #0, self.nbins - 1)
+                                      temp_first, temp_last)
+                                      #first_bin, last_bin)
+        if smear_output < 0:
+            msg = "_BaseSmearer: could not smear, code = %g" % smear_output
+            raise RuntimeError, msg
+
+        temp_first = first_bin + self.nbins_low
+        temp_last = self.nbins - self.nbins_high
+        out = iq_out[temp_first: temp_last]
+
+        return out
+    
+    def _initialize_smearer(self):
+        """
+        """
+        return NotImplemented
+            
+    
+    def _get_unextrapolated_bin(self, first_bin = 0, last_bin = 0):
+        """
+        Get unextrapolated first bin and the last bin
+        
+        : param first_bin: extrapolated first_bin
+        : param last_bin: extrapolated last_bin
+        
+        : return fist_bin, last_bin: unextrapolated first and last bin
+        """
+        #  For first bin
+        if first_bin <= self.nbins_low:
+            first_bin = 0
+        else:
+            first_bin = first_bin - self.nbins_low
+        # For last bin
+        if last_bin >= (self.nbins - self.nbins_high):
+            last_bin  = self.nbins - (self.nbins_high + self.nbins_low + 1)
+        elif last_bin >= self.nbins_low:
+            last_bin = last_bin - self.nbins_low
+        else:
+            last_bin = 0
+        return first_bin, last_bin
+    
+    def _get_extrapolated_bin(self, first_bin = 0, last_bin = 0):
+        """
+        Get extrapolated first bin and the last bin
+        
+        : param first_bin: unextrapolated first_bin
+        : param last_bin: unextrapolated last_bin
+        
+        : return first_bin, last_bin: extrapolated first and last bin
+        """
+        #  For the first bin
+        # In the case that needs low extrapolation data
+        first_bin = 0
+        # For last bin
+        if last_bin >= self.nbins - (self.nbins_high + self.nbins_low + 1):
+            # In the case that needs higher q extrapolation data 
+            last_bin = self.nbins - 1
+        else:
+            # In the case that doesn't need higher q extrapolation data 
+             last_bin += self.nbins_low
+
+        return first_bin, last_bin
+        
+class _SlitSmearer(_BaseSmearer):
+    """
+    Slit smearing for I(q) array
+    """
+    
+    def __init__(self, nbins=None, width=None, height=None, min=None, max=None):
+        """
+        Initialization
+            
+        :param iq: I(q) array [cm-1]
+        :param width: slit width [A-1]
+        :param height: slit height [A-1]
+        :param min: Q_min [A-1]
+        :param max: Q_max [A-1]
+        
+        """
+        _BaseSmearer.__init__(self)
+        ## Slit width in Q units
+        self.width  = width
+        ## Slit height in Q units
+        self.height = height
+        ## Q_min (Min Q-value for I(q))
+        self.min    = min
+        ## Q_max (Max Q_value for I(q))
+        self.max    = max
+        ## Number of Q bins 
+        self.nbins  = nbins
+        ## Number of points used in the smearing computation
+        self.npts   = 3000
+        ## Smearing matrix
+        self._weights = None
+        self.qvalues  = None
+        
+    def _initialize_smearer(self):
+        """
+        Initialize the C++ smearer object.
+        This method HAS to be called before smearing
+        """
+        #self._smearer = smearer.new_slit_smearer(self.width,
+        # self.height, self.min, self.max, self.nbins)
+        self._smearer = smearer.new_slit_smearer_with_q(self.width, 
+                                                    self.height, self.qvalues)
+        self._init_complete = True
+
+    def get_unsmeared_range(self, q_min, q_max):
+        """
+        Determine the range needed in unsmeared-Q to cover
+        the smeared Q range
+        """
+        # Range used for input to smearing
+        _qmin_unsmeared = q_min
+        _qmax_unsmeared = q_max 
+        try:
+            _qmin_unsmeared = self.min
+            _qmax_unsmeared = self.max
+        except:
+            logging.error("_SlitSmearer.get_bin_range: %s" % sys.exc_value)
+        return _qmin_unsmeared, _qmax_unsmeared
+
+class SlitSmearer(_SlitSmearer):
+    """
+    Adaptor for slit smearing class and SANS data
+    """
+    def __init__(self, data1D, model = None):
+        """
+        Assumption: equally spaced bins of increasing q-values.
+        
+        :param data1D: data used to set the smearing parameters
+        """
+        # Initialization from parent class
+        super(SlitSmearer, self).__init__()
+        
+        ## Slit width
+        self.width = 0
+        self.nbins_low = 0
+        self.nbins_high = 0
+        self.model = model
+        if data1D.dxw is not None and len(data1D.dxw) == len(data1D.x):
+            self.width = data1D.dxw[0]
+            # Sanity check
+            for value in data1D.dxw:
+                if value != self.width:
+                    msg = "Slit smearing parameters must "
+                    msg += " be the same for all data"
+                    raise RuntimeError, msg
+        ## Slit height
+        self.height = 0
+        if data1D.dxl is not None and len(data1D.dxl) == len(data1D.x):
+            self.height = data1D.dxl[0]
+            # Sanity check
+            for value in data1D.dxl:
+                if value != self.height:
+                    msg = "Slit smearing parameters must be"
+                    msg += " the same for all data"
+                    raise RuntimeError, msg
+        # If a model is given, get the q extrapolation
+        if self.model == None:
+            data1d_x = data1D.x
+        else:
+            # Take larger sigma
+            if self.height > self.width:
+                # The denominator (2.0) covers all the possible w^2 + h^2 range
+                sigma_in = data1D.dxl / 2.0
+            elif self.width > 0:
+                sigma_in = data1D.dxw / 2.0
+            else:
+                sigma_in = []
+
+            self.nbins_low, self.nbins_high, _, data1d_x = \
+                                get_qextrapolate(sigma_in, data1D.x)
+
+        ## Number of Q bins
+        self.nbins = len(data1d_x)
+        ## Minimum Q 
+        self.min = min(data1d_x)
+        ## Maximum
+        self.max = max(data1d_x)
+        ## Q-values
+        self.qvalues = data1d_x
+        
+
+class _QSmearer(_BaseSmearer):
+    """
+    Perform Gaussian Q smearing
+    """
+        
+    def __init__(self, nbins=None, width=None, min=None, max=None):
+        """
+        Initialization
+        
+        :param nbins: number of Q bins
+        :param width: array standard deviation in Q [A-1]
+        :param min: Q_min [A-1]
+        :param max: Q_max [A-1]
+        """
+        _BaseSmearer.__init__(self)
+        ## Standard deviation in Q [A-1]
+        self.width = width
+        ## Q_min (Min Q-value for I(q))
+        self.min = min
+        ## Q_max (Max Q_value for I(q))
+        self.max = max
+        ## Number of Q bins 
+        self.nbins = nbins
+        ## Smearing matrix
+        self._weights = None
+        self.qvalues  = None
+        
+    def _initialize_smearer(self):
+        """
+        Initialize the C++ smearer object.
+        This method HAS to be called before smearing
+        """
+        #self._smearer = smearer.new_q_smearer(numpy.asarray(self.width),
+        # self.min, self.max, self.nbins)
+        self._smearer = smearer.new_q_smearer_with_q(numpy.asarray(self.width),
+                                                      self.qvalues)
+        self._init_complete = True
+        
+    def get_unsmeared_range(self, q_min, q_max):
+        """
+        Determine the range needed in unsmeared-Q to cover
+        the smeared Q range
+        Take 3 sigmas as the offset between smeared and unsmeared space
+        """
+        # Range used for input to smearing
+        _qmin_unsmeared = q_min
+        _qmax_unsmeared = q_max 
+        try:
+            offset = 3.0 * max(self.width)
+            _qmin_unsmeared = self.min#max([self.min, q_min - offset])
+            _qmax_unsmeared = self.max#min([self.max, q_max + offset])
+        except:
+            logging.error("_QSmearer.get_bin_range: %s" % sys.exc_value)
+        return _qmin_unsmeared, _qmax_unsmeared
+        
+    
+class QSmearer(_QSmearer):
+    """
+    Adaptor for Gaussian Q smearing class and SANS data
+    """
+    def __init__(self, data1D, model = None):
+        """
+        Assumption: equally spaced bins of increasing q-values.
+        
+        :param data1D: data used to set the smearing parameters
+        """
+        # Initialization from parent class
+        super(QSmearer, self).__init__()
+        data1d_x = []
+        self.nbins_low = 0
+        self.nbins_high = 0
+        self.model = model
+        ## Resolution
+        #self.width = numpy.zeros(len(data1D.x))
+        if data1D.dx is not None and len(data1D.dx) == len(data1D.x):
+            self.width = data1D.dx
+        
+        if self.model == None:
+            data1d_x = data1D.x
+        else:
+            self.nbins_low, self.nbins_high, self.width, data1d_x = \
+                                get_qextrapolate(self.width, data1D.x)
+
+        ## Number of Q bins
+        self.nbins = len(data1d_x)
+        ## Minimum Q 
+        self.min = min(data1d_x)
+        ## Maximum
+        self.max = max(data1d_x)
+        ## Q-values
+        self.qvalues = data1d_x
+
+        
+def get_qextrapolate(width, data_x):
+    """
+    Make fake data_x points extrapolated outside of the data_x points
+    
+    : param width: array of std of q resolution
+    : param Data1D.x: Data1D.x array
+    
+    : return new_width, data_x_ext: extrapolated width array and x array
+    
+    : assumption1: data_x is ordered from lower q to higher q
+    : assumption2: len(data) = len(width)
+    : assumption3: the distance between the data points is more compact 
+            than the size of width 
+    : Todo1: Make sure that the assumptions are correct for Data1D
+    : Todo2: This fixes the edge problem in Qsmearer but still needs to make 
+            smearer interface 
+    """
+    # Length of the width
+    length = len(width)
+    width_low = math.fabs(width[0])
+    width_high = math.fabs(width[length -1])
+    
+    # Compare width(dQ) to the data bin size and take smaller one as the bin 
+    # size of the extrapolation; this will correct some weird behavior 
+    # at the edge: This method was out (commented) 
+    # because it becomes very expansive when
+    # bin size is very small comparing to the width.
+    # Now on, we will just give the bin size of the extrapolated points 
+    # based on the width.
+    # Find bin sizes
+    #bin_size_low = math.fabs(data_x[1] - data_x[0])
+    #bin_size_high = math.fabs(data_x[length - 1] - data_x[length - 2])
+    # Let's set the bin size 1/3 of the width(sigma), it is good as long as
+    # the scattering is monotonous.
+    #if width_low < (bin_size_low):
+    bin_size_low = width_low / 10.0
+    #if width_high < (bin_size_high):
+    bin_size_high = width_high / 10.0
+        
+    # Number of q points required below the 1st data point in order to extend
+    # them 3 times of the width (std)
+    nbins_low = math.ceil(3.0 * width_low / bin_size_low)
+    # Number of q points required above the last data point
+    nbins_high = math.ceil(3.0 * width_high / (bin_size_high))
+    # Make null q points        
+    extra_low = numpy.zeros(nbins_low)
+    extra_high = numpy.zeros(nbins_high)
+    # Give extrapolated values
+    ind = 0
+    qvalue = data_x[0] - bin_size_low
+    #if qvalue > 0:
+    while(ind < nbins_low):
+        extra_low[nbins_low - (ind + 1)] = qvalue
+        qvalue -= bin_size_low
+        ind += 1
+        #if qvalue <= 0:
+        #    break
+    # Redefine nbins_low
+    nbins_low = ind
+    # Reset ind for another extrapolation
+    ind = 0
+    qvalue = data_x[length -1] + bin_size_high
+    while(ind < nbins_high):
+        extra_high[ind] = qvalue
+        qvalue += bin_size_high
+        ind += 1
+    # Make a new qx array
+    if nbins_low > 0:  
+        data_x_ext = numpy.append(extra_low, data_x)
+    else:
+        data_x_ext = data_x
+    data_x_ext = numpy.append(data_x_ext, extra_high)
+    
+    # Redefine extra_low and high based on corrected nbins  
+    # And note that it is not necessary for extra_width to be a non-zero 
+    if nbins_low > 0:     
+        extra_low = numpy.zeros(nbins_low)
+    extra_high = numpy.zeros(nbins_high) 
+    # Make new width array
+    new_width = numpy.append(extra_low, width)
+    new_width = numpy.append(new_width, extra_high)
+    
+    # nbins corrections due to the negative q value
+    nbins_low = nbins_low - len(data_x_ext[data_x_ext<=0])
+    return  nbins_low, nbins_high, \
+             new_width[data_x_ext>0], data_x_ext[data_x_ext>0]
+    
+if __name__ == '__main__':
+    x = 0.001 * numpy.arange(1, 11)
+    y = 12.0 - numpy.arange(1, 11)
+    print x
+    #for i in range(10): print i, 0.001 + i*0.008/9.0 
+    #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) )) 
+    s = _SlitSmearer(nbins=10, width=0.0, height=0.005, min=0.001, max=0.010)
+    #s = _QSmearer(nbins=10, width=0.001, min=0.001, max=0.010)
+    s._compute_matrix()
+
+    sy = s(y)
+    print sy
+    
+    if True:
+        for i in range(10):
+            print x[i], y[i], sy[i]
+            #print q, ' : ', s.weight(q), s._compute_iq(q) 
+            #print q, ' : ', s(q), s._compute_iq(q) 
+            #s._compute_iq(q) 
+
+
+
+
Index: sansmodels/src/sans/models/smearing_2d.py
===================================================================
--- sansmodels/src/sans/models/smearing_2d.py	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
+++ sansmodels/src/sans/models/smearing_2d.py	(revision 642b259b9abb82689ce5f576e23b9b02fdb74600)
@@ -0,0 +1,255 @@
+#####################################################################
+#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 2008, University of Tennessee
+######################################################################
+
+## TODO: Need test,and check Gaussian averaging
+import numpy, math,time
+## Singular point
+SIGMA_ZERO = 1.0e-010
+## Limit of how many sigmas to be covered for the Gaussian smearing
+# default: 2.5 to cover 98.7% of Gaussian
+LIMIT = 2.5
+## 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}   
+
+class Smearer2D:
+    """
+    Gaussian Q smearing class for SANS 2d data
+    """
+     
+    def __init__(self, data=None, model=None, index=None, 
+                 limit=LIMIT, accuracy='Low'):
+        """
+        Assumption: equally spaced bins in dq_r, dq_phi space.
+        
+        :param data: 2d data used to set the smearing parameters
+        :param model: model function
+        :param index: 1d array with len(data) to define the range of the calculation: elements are given as True or False
+        :param nr: number of bins in dq_r-axis
+        :param nphi: number of bins in dq_phi-axis
+        """
+        ## data
+        self.data = data
+        ## model
+        self.model = model
+        ## Accuracy: Higher stands for more sampling points in both directions of r and phi.
+        self.accuracy = accuracy
+        ## number of bins in r axis for over-sampling 
+        self.nr = R_BIN
+        ## number of bins in phi axis for over-sampling 
+        self.nphi = PHI_BIN
+        ## maximum nsigmas
+        self.limit = limit
+        self.index = index
+        self.smearer = True
+        
+        
+    def get_data(self):   
+        """
+        get qx_data, qy_data, dqx_data,dqy_data,and calculate phi_data=arctan(qx_data/qy_data)
+        """
+        if self.data == None or self.data.__class__.__name__ == 'Data1D':
+            return None
+        if self.data.dqx_data == None or self.data.dqy_data == None:
+            return None
+        self.qx_data = self.data.qx_data[self.index]
+        self.qy_data = self.data.qy_data[self.index]
+        self.dqx_data = self.data.dqx_data[self.index]
+        self.dqy_data = self.data.dqy_data[self.index]
+        self.phi_data = numpy.arctan(self.qx_data/self.qy_data)
+        ## Remove singular points if exists
+        self.dqx_data[self.dqx_data<SIGMA_ZERO]=SIGMA_ZERO
+        self.dqy_data[self.dqy_data<SIGMA_ZERO]=SIGMA_ZERO
+        return True
+    
+    def set_accuracy(self, accuracy='Low'):   
+        """
+        Set accuracy.
+        
+        :param accuracy:  string
+        """
+        self.accuracy = accuracy
+
+    def set_smearer(self, smearer=True):
+        """
+        Set whether or not smearer will be used
+        
+        :param smearer: smear object
+        
+        """
+        self.smearer = smearer
+        
+    def set_data(self, data=None):   
+        """
+        Set data.
+        
+        :param data: DataLoader.Data_info type
+        """
+        self.data = data
+  
+            
+    def set_model(self, model=None):   
+        """
+        Set model.
+        
+        :param model: sans.models instance
+        """
+        self.model = model  
+           
+    def set_index(self, index=None):   
+        """
+        Set index.
+        
+        :param index: 1d arrays
+        """
+        self.index = index       
+    
+    def get_value(self):
+        """
+        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:
+            return valid
+        if self.smearer == False:
+            return self.model.evalDistribution([self.qx_data,self.qy_data]) 
+        st = time.time()
+        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
+        # Total number of bins = # of bins 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)+numpy.arctan(self.qy_data*self.dqx_data/self.qx_data/self.dqy_data).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)))
+        weight_res = weight_res.repeat(nphi).reshape(nr,nphi).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()
+
+        
+        ## Over-sampled qx_data and qy_data.
+        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
+    
+if __name__ == '__main__':
+    ## Test w/ 2D linear function
+    x = 0.001*numpy.arange(1,11)
+    dx = numpy.ones(len(x))*0.001
+    y = 0.001*numpy.arange(1,11)
+    dy = numpy.ones(len(x))*0.001
+    z = numpy.ones(10)
+    dz = numpy.sqrt(z)
+    
+    from DataLoader import Data2D
+    #for i in range(10): print i, 0.001 + i*0.008/9.0 
+    #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) )) 
+    out = Data2D()
+    out.data = z
+    out.qx_data = x
+    out.qy_data = y
+    out.dqx_data = dx
+    out.dqy_data = dy
+    index = numpy.ones(len(x), dtype = bool)
+    out.mask = index
+    from sans.models.LineModel import LineModel
+    model = LineModel()
+    model.setParam("A", 0)
+
+    smear = Smearer2D(out,model,index)
+    #smear.set_accuracy('Xhigh')
+    value = smear.get_value()
+    ## All data are ones, so the smeared should also be ones.
+    print "Data length =",len(value)
+    print " 2D linear function, I = 0 + 1*qx*qy"
+    print " Gaussian weighted averaging on a 2D linear function will provides the results same as without the averaging."
+    print "qx_data", "qy_data", "I_nonsmear", "I_smeared"
+    for ind in range(len(value)):
+        print x[ind],y[ind],model.evalDistribution([x,y])[ind], value[ind]
+  
+"""    
+if __name__ == '__main__':
+    ## Another Test w/ constant function
+    x = 0.001*numpy.arange(1,11)
+    dx = numpy.ones(len(x))*0.001
+    y = 0.001*numpy.arange(1,11)
+    dy = numpy.ones(len(x))*0.001
+    z = numpy.ones(10)
+    dz = numpy.sqrt(z)
+    
+    from DataLoader import Data2D
+    #for i in range(10): print i, 0.001 + i*0.008/9.0 
+    #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) )) 
+    out = Data2D()
+    out.data = z
+    out.qx_data = x
+    out.qy_data = y
+    out.dqx_data = dx
+    out.dqy_data = dy
+    index = numpy.ones(len(x), dtype = bool)
+    out.mask = index
+    from sans.models.Constant import Constant
+    model = Constant()
+
+    value = Smearer2D(out,model,index).get_value()
+    ## All data are ones, so the smeared values should also be ones.
+    print "Data length =",len(value), ", Data=",value
+"""    
