Index: DataLoader/extensions/smearer.cpp
===================================================================
--- DataLoader/extensions/smearer.cpp	(revision a74d650ab9d64fbf12a4d8a32772eb925dfaa913)
+++ DataLoader/extensions/smearer.cpp	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
@@ -133,64 +133,134 @@
 	// 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++) {
-		double q, q_min, q_max;
+		// 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);
-
-		// For each q-value, compute the weight of each other q-bin
-		// in the I(q) array
-		int npts_h = height>0 ? npts : 1;
-		int npts_w = width>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)
-		// Never down-grade npts_h. That will give incorrect slit smearing...
-		if(npts_h>1 && npts_w>1){
-			npts_h = npts;//(int)ceil(sqrt((double)npts));
-			// In general width is much smaller than height, so smaller npt_w
-			// should work fine.
-			// Todo: It is still very expansive in time. Think about better way.
-			npts_w = (int)ceil(npts_h / 100);
-		}
-
-		double shift_h, shift_w;
-		for(int k=0; k<npts_h; k++){
-			if(npts_h==1){
-				shift_h = 0;
-			} else {
-				shift_h = height/((double)npts_h-1.0) * (double)k;
-			}
-			for(int j=0; j<npts_w; j++){
-				if(npts_w==1){
-					shift_w = 0;
-				} else {
-					shift_w = width/((double)npts_w-1.0) * (double)j;
+		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);
+			// 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;
+			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;
 				}
-				double q_shifted = sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) );
-
-				// Find in which bin this shifted value belongs
-				int q_i=nbins;
-				if (even_binning) {
-					// This is kept for backward compatibility since the binning
-					// was originally defined differently for even bins.
-					q_i = (int)(floor( (q_shifted-qmin) /((qmax-qmin)/((double)nbins -1.0)) ));
-				} else {
-					for(int t=0; t<nbins; t++) {
-						double q_t, q_high, q_low;
-						get_bin_range(t, &q_t, &q_low, &q_high);
-						if(q_shifted>=q_low && q_shifted<q_high) {
-							q_i = t;
-							break;
+			}
+			//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 ;
 						}
 					}
 				}
-
-				// Skip the entries outside our I(q) range
-				//TODO: be careful with edge effect
-				if(q_i<nbins)
-					(*weights)[i*nbins+q_i]++;
 			}
 		}
@@ -277,8 +347,8 @@
 
 		for(int i=first_bin; i<=last_bin; i++){
-			// Skip if weight is less than 1e-04(this value is much smaller than
+			// 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-004){
+			if ((*weights)[q_i*nbins+i] < 1.0e-003){
 				continue;
 			}
@@ -288,5 +358,5 @@
 
 		// Normalize counts
-		iq_out[q_i] = (counts>0.0) ? sum/counts : 0;
+		iq_out[q_i] = (counts>0.0) ? sum/counts : 0.0;
 	}
 }
Index: DataLoader/extensions/smearer.hh
===================================================================
--- DataLoader/extensions/smearer.hh	(revision 2ca00c4c8465dbd97c9ce82a83b84c24dd55d1a6)
+++ DataLoader/extensions/smearer.hh	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
@@ -61,5 +61,5 @@
 protected:
     // Number of points used in the smearing computation
-    static const int npts   = 3000;
+    static const int npts   = 500;
 
 public:
Index: DataLoader/qsmearing.py
===================================================================
--- DataLoader/qsmearing.py	(revision 2ca00c4c8465dbd97c9ce82a83b84c24dd55d1a6)
+++ DataLoader/qsmearing.py	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
@@ -157,4 +157,6 @@
         Perform smearing
         """
+        import time
+        st = time.time()
         # If this is the first time we call for smearing,
         # initialize the C++ smearer object first
@@ -174,5 +176,5 @@
             temp_first, temp_last = self._get_extrapolated_bin( \
                                                         first_bin, last_bin)
-            if self.nbins_low > 1:
+            if self.nbins_low > 0:
                 iq_in_low = self.model.evalDistribution( \
                                     numpy.fabs(self.qvalues[0:self.nbins_low]))
@@ -213,5 +215,5 @@
         temp_last = self.nbins - self.nbins_high
         out = iq_out[temp_first: temp_last]
-
+        print "time =", time.time() - st
         return out
     
@@ -549,7 +551,7 @@
     # Make a new qx array
     if nbins_low > 0:  
-    	data_x_ext = numpy.append(extra_low, data_x)
+        data_x_ext = numpy.append(extra_low, data_x)
     else:
-    	data_x_ext = data_x
+        data_x_ext = data_x
     data_x_ext = numpy.append(data_x_ext, extra_high)
     
Index: DataLoader/test/slit_1000A_sphere_sm_h.txt
===================================================================
--- DataLoader/test/slit_1000A_sphere_sm_h.txt	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
+++ DataLoader/test/slit_1000A_sphere_sm_h.txt	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
@@ -0,0 +1,76 @@
+<X>   <Y>
+3.0525e-005  272927
+3.33e-005  271895
+3.885e-005  269570
+4.44e-005  266907
+4.995e-005  263918
+5.55e-005  260610
+6.105e-005  256996
+6.66e-005  253088
+7.215e-005  248898
+7.77e-005  244441
+8.325e-005  239731
+8.88e-005  234783
+9.435e-005  229615
+9.99e-005  224244
+0.00010545  218691
+0.000111  213007
+0.0001221  201099
+0.0001332  188720
+0.0001443  176014
+0.0001554  163125
+0.0001665  150194
+0.0001776  137357
+0.0001887  124740
+0.0001998  112461
+0.0002109  100628
+0.000222  89334.1
+0.0002331  78661.2
+0.0002442  68678.3
+0.0002553  59444.6
+0.0002664  51021.4
+0.0002775  43572.8
+0.00030525  27564.1
+0.000333  16260.5
+0.00036075  9066.14
+0.0003885  5134.76
+0.00041625  3520.4
+0.000444  3314.21
+0.00047175  3749.34
+0.0004995  4263.86
+0.00052725  4519.32
+0.000555  4377.13
+0.0006105  3267.94
+0.000666  1795.2
+0.0007215  914.248
+0.000777  792.066
+0.0008325  981.327
+0.000888  1011.67
+0.0009435  772.856
+0.000999  462.554
+0.0010545  296.986
+0.00111  304.392
+0.001221  372.281
+0.001332  200.069
+0.001443  155.226
+0.001554  174.773
+0.001665  95.527
+0.001776  91.3966
+0.001887  93.2306
+0.001998  53.2604
+0.002109  58.6666
+0.00222  53.9608
+0.002331  33.3167
+0.002442  39.6982
+0.002553  33.0803
+0.002664  22.8073
+0.002775  28.1541
+0.0030525  14.841
+0.00333  9.58852
+0.0036075  9.25752
+0.003885  10.6765
+0.0041625  11.6735
+0.00444  11.3733
+0.0047175  9.95799
+0.004995  8.1725
+0.0052725  6.8853
Index: DataLoader/test/slit_1000A_sphere_sm_w_0_0001.txt
===================================================================
--- DataLoader/test/slit_1000A_sphere_sm_w_0_0001.txt	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
+++ DataLoader/test/slit_1000A_sphere_sm_w_0_0001.txt	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
@@ -0,0 +1,76 @@
+<X>   <Y>
+3.0525e-005  1.58586e+008
+3.33e-005  1.55688e+008
+3.885e-005  1.55688e+008
+4.44e-005  1.52645e+008
+4.995e-005  1.52645e+008
+5.55e-005  1.49479e+008
+6.105e-005  1.49479e+008
+6.66e-005  1.46208e+008
+7.215e-005  1.46208e+008
+7.77e-005  1.42852e+008
+8.325e-005  1.42852e+008
+8.88e-005  1.3943e+008
+9.435e-005  1.3943e+008
+9.99e-005  1.35962e+008
+0.00010545  1.35361e+008
+0.000111  1.29681e+008
+0.0001221  1.23716e+008
+0.0001332  1.17819e+008
+0.0001443  1.11547e+008
+0.0001554  1.04755e+008
+0.0001665  9.78708e+007
+0.0001776  8.86188e+007
+0.0001887  8.50559e+007
+0.0001998  8.14237e+007
+0.0002109  7.00377e+007
+0.000222  6.72838e+007
+0.0002331  5.69979e+007
+0.0002442  5.35122e+007
+0.0002553  5.00696e+007
+0.0002664  4.13422e+007
+0.0002775  3.81951e+007
+0.00030525  2.55485e+007
+0.000333  1.78886e+007
+0.00036075  1.0397e+007
+0.0003885  5.5246e+006
+0.00041625  2.86391e+006
+0.000444  1.42479e+006
+0.00047175  841250
+0.0004995  632587
+0.00052725  789548
+0.000555  909313
+0.0006105  988789
+0.000666  641836
+0.0007215  251620
+0.000777  80779.6
+0.0008325  105819
+0.000888  172064
+0.0009435  170644
+0.000999  104952
+0.0010545  34613.8
+0.00111  8640.37
+0.001221  66520.7
+0.001332  23138.3
+0.001443  4608.17
+0.001554  27002.8
+0.001665  6005.28
+0.001776  4189.92
+0.001887  12470.6
+0.001998  1443.87
+0.002109  3468.11
+0.00222  6127.33
+0.002331  231.738
+0.002442  2735.44
+0.002553  3069.6
+0.002664  2.14309
+0.002775  2087.64
+0.0030525  781.373
+0.00333  146.203
+0.0036075  0.623872
+0.003885  102.714
+0.0041625  252.086
+0.00444  334.016
+0.0047175  318.556
+0.004995  233.168
+0.0052725  148.507
Index: DataLoader/test/slit_1000A_sphere_sm_w_0_0002.txt
===================================================================
--- DataLoader/test/slit_1000A_sphere_sm_w_0_0002.txt	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
+++ DataLoader/test/slit_1000A_sphere_sm_w_0_0002.txt	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
@@ -0,0 +1,76 @@
+<X>   <Y>
+3.0525e-005  213022
+3.33e-005  217281
+3.885e-005  215212
+4.44e-005  213920
+4.995e-005  212464
+5.55e-005  210809
+6.105e-005  209013
+6.66e-005  207049
+7.215e-005  204968
+7.77e-005  202733
+8.325e-005  200231
+8.88e-005  197676
+9.435e-005  195074
+9.99e-005  192314
+0.00010545  189393
+0.000111  196598
+0.0001221  188633
+0.0001332  181984
+0.0001443  174955
+0.0001554  167657
+0.0001665  160010
+0.0001776  152171
+0.0001887  144150
+0.0001998  135955
+0.0002109  127651
+0.000222  119142
+0.0002331  110131
+0.0002442  101201
+0.0002553  92502.5
+0.0002664  84098.4
+0.0002775  82797.3
+0.00030525  61108.1
+0.000333  44383
+0.00036075  31119
+0.0003885  21003.7
+0.00041625  13791.8
+0.000444  8958.24
+0.00047175  5995.06
+0.0004995  4537.52
+0.00052725  3924.02
+0.000555  4049
+0.0006105  3462.35
+0.000666  2543.48
+0.0007215  1557.49
+0.000777  994.03
+0.0008325  870.131
+0.000888  893.879
+0.0009435  820.013
+0.000999  623.98
+0.0010545  412.198
+0.00111  316.687
+0.001221  360.124
+0.001332  236.643
+0.001443  145.266
+0.001554  174.021
+0.001665  110.306
+0.001776  83.4758
+0.001887  95.9231
+0.001998  59.2528
+0.002109  53.4788
+0.00222  57.1684
+0.002331  35.3826
+0.002442  36.6378
+0.002553  35.9109
+0.002664  21.7072
+0.002775  28.337
+0.0030525  14.8856
+0.00333  9.61047
+0.0036075  9.32256
+0.003885  10.791
+0.0041625  11.81
+0.00444  11.4979
+0.0047175  10.0494
+0.004995  8.2325
+0.0052725  6.88521
Index: DataLoader/test/utest_smearing.py
===================================================================
--- DataLoader/test/utest_smearing.py	(revision f867cd99b41e1e3c36ce344d9791e93e37fe6bde)
+++ DataLoader/test/utest_smearing.py	(revision cd2ced800865b1cac7a3bcda776972609f33927f)
@@ -9,5 +9,5 @@
 from DataLoader.data_info import Data1D, Data2D
 from DataLoader.qsmearing import SlitSmearer, QSmearer, smear_selection
- 
+from sans.models.SphereModel import SphereModel
 import os.path
 
@@ -54,7 +54,15 @@
         input = 12.0-numpy.arange(1,11)
         output = s(input)
-        # The following commented line was the correct output for even bins [see smearer.cpp for details] 
-        #answer = [ 9.666,  9.056,  8.329,  7.494,  6.642,  5.721,  4.774,  3.824,  2.871, 2.   ]
-        answer = [ 9.2302,  8.6806,  7.9533,  7.1673,  6.2889,  5.4,     4.5028,  3.5744,  2.6083, 2.    ]
+        # The following commented line was the correct output for even bins 
+        # [see smearer.cpp for details] 
+        #answer = [ 9.666,  9.056,  8.329,  7.494,  6.642,  5.721,  4.774, \
+        #  3.824,  2.871, 2.   ]
+        # The following answer was from numerical weighting algorithm.
+        #answer = [ 9.2302,  8.6806,  7.9533,  7.1673,  6.2889,  5.4,   \
+        #  4.5028,  3.5744,  2.6083, 2.    ]
+        # For the new analytical algorithm, the small difference between 
+        #these two could be from the first edge of the q bin size. 
+        answer = [ 9.0618,  8.64018,  8.11868,  7.13916,  6.15285,  5.55556,  \
+                     4.55842,  3.56061,  2.56235, 2.    ]
         for i in range(len(input)):
             self.assertAlmostEqual(answer[i], output[i], 2)
@@ -89,7 +97,82 @@
                   5.00093415,   4.01898292,   3.15008701,   2.55214921]
         for i in range(len(input)):
-            self.assertAlmostEqual(answer[i], output[i], 4)
-      
-
+            self.assertAlmostEqual(answer[i], output[i], 2)
+            
+class smear_slit_h_w_tests(unittest.TestCase):
+    
+    def setUp(self):
+        self.data = Loader().load("1000A_sphere_sm.xml")
+        self.model = SphereModel()
+        # The answer could be improved by developing better algorithm.
+        self.answer1 = Loader().load("slit_1000A_sphere_sm_w_0_0002.txt")
+        self.answer2 = Loader().load("slit_1000A_sphere_sm_h.txt")
+        self.answer3 = Loader().load("slit_1000A_sphere_sm_w_0_0001.txt")
+        # Get inputs
+        self.model.params['scale'] = 0.05
+        self.model.params['background'] = 0.01
+        self.model.params['radius'] = 10000.0
+        self.model.params['sldSolv'] = 6.3e-006
+        self.model.params['sldSph'] = 3.4e-006
+        
+    def test_slit_h_w(self):
+        """
+            Test identity slit smearing w/ h=0.117 w = 0.002
+        """
+        # Set params and dQl
+        data = self.data
+        data.dxw = 0.0002 * numpy.ones(len(self.data.x))
+        data.dxl = 0.117 * numpy.ones(len(self.data.x))
+        # Create smearer for our data
+        s = SlitSmearer(data, self.model)
+        # Get smear
+        input  = self.model.evalDistribution(data.x)
+        output = s(input)
+        # Get pre-calculated values
+        answer = self.answer1.y
+        # Now compare
+        for i in range(len(input)):
+            self.assertAlmostEqual(answer[i], output[i], 0)
+     
+    def test_slit_h(self):
+        """
+            Test identity slit smearing w/ h=0.117 w = 0.0
+        """
+        # Set params and dQl
+        data = self.data
+        data.dxw = 0.0 * numpy.ones(len(self.data.x))
+        data.dxl = 0.117 * numpy.ones(len(self.data.x))
+        # Create smearer for our data
+        s = SlitSmearer(data, self.model)
+        # Get smear
+        input  = self.model.evalDistribution(data.x)
+        output = s(input)
+        # Get pre-calculated values
+        answer = self.answer2.y
+        # Now compare
+        for i in range(len(input)):
+            self.assertAlmostEqual(answer[i], output[i], 0)
+     
+    def test_slit_w(self):
+        """
+            Test identity slit smearing w/ h=0.0 w = 0.001
+        """
+        # Set params and dQl
+        data = self.data
+        data.dxw = 0.0001 * numpy.ones(len(self.data.x))
+        data.dxl = 0.0 * numpy.ones(len(self.data.x))
+        # Create smearer for our data
+        s = SlitSmearer(data, self.model)
+        # Get smear
+        input  = self.model.evalDistribution(data.x)
+        output = s(input)
+        # Get pre-calculated values
+        answer = self.answer3.y
+        # Now compare
+        for i in range(len(input)):
+             if i <= 40:
+                 self.assertAlmostEqual(answer[i], output[i], -3)    
+             else:
+                self.assertAlmostEqual(answer[i], output[i], 0)       
+                
 if __name__ == '__main__':
     unittest.main()
