Index: sansmodels/src/sans/models/BaseComponent.py
===================================================================
--- sansmodels/src/sans/models/BaseComponent.py	(revision 26e4a247fb629283560bcebfc702766d91e89b16)
+++ sansmodels/src/sans/models/BaseComponent.py	(revision ecc58e7254a3e176f1f68af184d95d2ec777177a)
@@ -60,20 +60,68 @@
         """
             Evaluate a distribution of q-values.
-            A list of either scale q-values or [qx,qy] doublets
-            are assumed. The allowed data types for the doublets are 
-            the same as for run() and runXY().
-            
-            Since the difference between scale and vector q-values
-            are dealt with in runXY(), we pass along the values
-            as-is.
-            
-            @param qdist: list of scalar q-values or list of [qx,qy] doublets 
-        """
-        q_array = numpy.asarray(qdist)
-        iq_array = numpy.zeros(len(q_array))
-        for i in xrange(len(q_array)):
-            iq_array[i] = self.runXY(q_array[i])
-            
-        return iq_array
+            
+            * For 1D, a numpy array is expected as input:
+            
+                evalDistribution(q)
+                
+            where q is a numpy array.
+            
+            
+            * For 2D, a list of numpy arrays are expected: [qx_prime,qy_prime].
+            
+            For the following matrix of q_values [qx,qy], we will have:
+                      +--------+--------+--------+
+            qy[2]     |        |        |        | 
+                      +--------+--------+--------+
+            qy[1]     |        |        |        |
+                      +--------+--------+--------+            
+            qy[0]     |        |        |        | 
+                      +--------+--------+--------+
+                        qx[0]    qx[1]    qx[2]
+
+
+            The q-values must be stored in numpy arrays.
+            The qxprime and qyprime are 2D arrays. From the two lists of q-values 
+            above, numpy arrays must be created the following way:
+
+            qx_prime = numpy.reshape(qx, [3,1])
+            qy_prime = numpy.reshape(qy, [1,3])
+            
+            The method is then called the following way:
+            
+            evalDistribution([qx_prime, qy_prime])
+            
+            @param qdist: ndarray of scalar q-values or list [qx,qy] where qx,qy are 2D ndarrays 
+        """
+        if qdist.__class__.__name__ == 'list':
+            # Check whether we have a list of ndarrays [qx,qy]
+            if len(qdist)!=2 or \
+                qdist[0].__class__.__name__ != 'ndarray' or \
+                qdist[1].__class__.__name__ != 'ndarray':
+                    raise RuntimeError, "evalDistribution expects a list of 2 ndarrays"
+                
+            # Extract qx and qy for code clarity
+            qx = qdist[0]
+            qy = qdist[1]
+            
+            # Create output array
+            iq_array = numpy.zeros([qx.shape[0], qy.shape[1]])
+             
+            for i in range(qx.shape[0]):
+                for j in range(qy.shape[1]):
+                    iq_array[i][j] = self.runXY([qx[i][0],qy[0][j]])
+            return iq_array
+                
+        elif qdist.__class__.__name__ == 'ndarray':
+                # We have a simple 1D distribution of q-values
+                iq_array = numpy.zeros(len(qdist))
+                for i in range(len(qdist)):
+                    iq_array[i] = self.runXY(qdist[i])
+                return iq_array
+            
+        else:
+            mesg = "evalDistribution is expecting an ndarray of scalar q-values"
+            mesg += " or a list [qx,qy] where qx,qy are 2D ndarrays."
+            raise RuntimeError, mesg
     
     def clone(self):
Index: sansmodels/src/sans/models/test/utest_evaldist.py
===================================================================
--- sansmodels/src/sans/models/test/utest_evaldist.py	(revision ecc58e7254a3e176f1f68af184d95d2ec777177a)
+++ sansmodels/src/sans/models/test/utest_evaldist.py	(revision ecc58e7254a3e176f1f68af184d95d2ec777177a)
@@ -0,0 +1,165 @@
+"""
+    Test for the BaseComponent.evalDistribution
+    See method documentation for more details
+"""
+import unittest
+import numpy, math
+from sans.models.SphereModel import SphereModel
+from sans.models.Cos import Cos
+
+
+class TestEvalPythonMethods(unittest.TestCase):
+    """ Testing evalDistribution for pure python models """
+
+    def setUp(self):
+        self.model= Cos()
+        
+    def test_scalar_methods(self):
+        """
+            Simple test comparing the run(), runXY() and
+            evalDistribution methods
+        """
+        q1 = self.model.run(0.001)
+        q2 = self.model.runXY(0.001)
+        qlist3 = numpy.asarray([0.001, 0.002])
+        q3 = self.model.evalDistribution(qlist3)
+        q4 = self.model.run(0.002)
+
+        self.assertEqual(q1, q2)
+        self.assertEqual(q1, q3[0])
+        self.assertEqual(q4, q3[1])
+        
+    def test_XY_methods(self):
+        """
+            Compare to the runXY() method for 2D models.
+            
+                      +--------+--------+--------+
+            qy=0.009  |        |        |        | 
+                      +--------+--------+--------+
+            qy-0.006  |        |        |        |
+                      +--------+--------+--------+            
+            qy=0.003  |        |        |        | 
+                      +--------+--------+--------+
+                      qx=0.001   0.002   0.003
+            
+        """
+        # These are the expected values for all bins
+        expected = numpy.zeros([3,3])
+        for i in range(3):
+            for j in range(3):
+                q_length = math.sqrt( (0.001*(i+1.0))*(0.001*(i+1.0)) + (0.003*(j+1.0))*(0.003*(j+1.0)) )
+                expected[i][j] = self.model.run(q_length)
+        
+        qx_values = [0.001, 0.002, 0.003]
+        qy_values = [0.003, 0.006, 0.009]
+        
+        qx = numpy.asarray(qx_values)
+        qy = numpy.asarray(qy_values)
+        
+        qx_prime = numpy.reshape(qx, [3,1])
+        qy_prime = numpy.reshape(qy, [1,3])
+        
+        iq = self.model.evalDistribution([qx_prime, qy_prime])
+
+        for i in range(3):
+            for j in range(3):
+                self.assertAlmostEquals(iq[i][j], expected[i][j])
+
+    def test_rectangle_methods(self):
+        """
+            Compare to the runXY() method for 2D models
+            with a non-square matrix.
+            TODO: Doesn't work for C models apparently
+        
+                      +--------+--------+--------+
+            qy-0.006  |        |        |        |
+                      +--------+--------+--------+            
+            qy=0.003  |        |        |        | 
+                      +--------+--------+--------+
+                      qx=0.001   0.002   0.003
+            
+        """
+        # These are the expected values for all bins
+        expected = numpy.zeros([3,3])
+        for i in range(3):
+            for j in range(2):
+                q_length = math.sqrt( (0.001*(i+1.0))*(0.001*(i+1.0)) + (0.003*(j+1.0))*(0.003*(j+1.0)) )
+                expected[i][j] = self.model.run(q_length)
+        
+        qx_values = [0.001, 0.002, 0.003]
+        qy_values = [0.003, 0.006]
+        
+        qx = numpy.asarray(qx_values)
+        qy = numpy.asarray(qy_values)
+        
+        qx_prime = numpy.reshape(qx, [3,1])
+        qy_prime = numpy.reshape(qy, [1,2])
+        
+        iq = self.model.evalDistribution([qx_prime, qy_prime])
+        
+        for i in range(3):
+            for j in range(2):
+                self.assertAlmostEquals(iq[i][j], expected[i][j])
+
+
+                
+class TestEvalMethods(unittest.TestCase):
+    """ Testing evalDistribution for C models """
+
+    def setUp(self):
+        self.model= SphereModel()
+        
+    def test_scalar_methods(self):
+        """
+            Simple test comparing the run(), runXY() and
+            evalDistribution methods
+        """
+        q1 = self.model.run(0.001)
+        q2 = self.model.runXY(0.001)
+        q4 = self.model.run(0.002)
+        qlist3 = numpy.asarray([0.001, 0.002])
+        q3 = self.model.evalDistribution(qlist3)
+
+        self.assertEqual(q1, q2)
+        self.assertEqual(q1, q3[0])
+        self.assertEqual(q4, q3[1])
+        
+    def test_XY_methods(self):
+        """
+            Compare to the runXY() method for 2D models.
+            
+                      +--------+--------+--------+
+            qy=0.009  |        |        |        | 
+                      +--------+--------+--------+
+            qy-0.006  |        |        |        |
+                      +--------+--------+--------+            
+            qy=0.003  |        |        |        | 
+                      +--------+--------+--------+
+                      qx=0.001   0.002   0.003
+            
+        """
+        # These are the expected values for all bins
+        expected = numpy.zeros([3,3])
+        for i in range(3):
+            for j in range(3):
+                q_length = math.sqrt( (0.001*(i+1.0))*(0.001*(i+1.0)) + (0.003*(j+1.0))*(0.003*(j+1.0)) )
+                expected[i][j] = self.model.run(q_length)
+        
+        qx_values = [0.001, 0.002, 0.003]
+        qy_values = [0.003, 0.006, 0.009]
+        
+        qx = numpy.asarray(qx_values)
+        qy = numpy.asarray(qy_values)
+        
+        qx_prime = numpy.reshape(qx, [3,1])
+        qy_prime = numpy.reshape(qy, [1,3])
+        
+        iq = self.model.evalDistribution([qx_prime, qy_prime])
+        
+        for i in range(3):
+            for j in range(3):
+                self.assertAlmostEquals(iq[i][j], expected[i][j])
+                
+
+if __name__ == '__main__':
+    unittest.main()
