Index: sansmodels/src/sans/models/c_extensions/c_models.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/c_models.c	(revision 2c4b289dd897e67655e7c208484ad09f086d1094)
+++ sansmodels/src/sans/models/c_extensions/c_models.c	(revision 812b901e8b8451f555bc0649e5ee190827ff9130)
@@ -6,5 +6,6 @@
  */
 #include <Python.h>
-
+#define PY_ARRAY_UNIQUE_SYMBOL PyArray_API_sans
+#include "arrayobject.h"
 /**
  * Define empty module
@@ -25,5 +26,6 @@
     m = Py_InitModule3("c_models", module_methods,
                        "C extension module for SANS scattering models.");
-
+	
+	import_array();
 	addCCylinderModel(m);
 	addCParallelepipedModel(m);
Index: sansmodels/src/sans/models/c_models/CCylinderModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CCylinderModel.cpp	(revision 0f5bc9ff2027432913fd8d2cb067675244acdaf6)
+++ sansmodels/src/sans/models/c_models/CCylinderModel.cpp	(revision 812b901e8b8451f555bc0649e5ee190827ff9130)
@@ -23,6 +23,10 @@
  */
  
+#define NO_IMPORT_ARRAY
+#define PY_ARRAY_UNIQUE_SYMBOL PyArray_API_sans
+
 extern "C" {
 #include <Python.h>
+#include <arrayobject.h>
 #include "structmember.h"
 #include <stdio.h>
@@ -142,6 +146,107 @@
     }
 }
-
-
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array q[] 
+ * @return: numpy array object 
+ */
+ 
+static PyObject *evaluateOneDim(CylinderModel* model, PyArrayObject *q){
+    PyArrayObject *result;
+   
+    // Check validity of array q , q must be of dimension 1, an array of double
+    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE)
+    {
+        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message);
+        return NULL;
+    }
+    result = (PyArrayObject *)PyArray_FromDims(q->nd, (int *)(q->dimensions), 
+										  PyArray_DOUBLE);
+	if (result == NULL) {
+        const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+		return NULL;
+	}
+	 for (int i = 0; i < q->dimensions[0]; i++){
+      double q_value  = *(double *)(q->data + i*q->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      *result_value =(*model)(q_value);
+	}
+    return PyArray_Return(result); 
+ }
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array  [q[],phi[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDim( CylinderModel* model, 
+                              PyArrayObject *q, PyArrayObject *phi)
+ {
+    PyArrayObject *result;
+    //check validity of input vectors
+    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE
+        || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE
+        || phi->dimensions[0] != q->dimensions[0]){
+        
+        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+	result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE);
+
+	if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	}
+	
+    for (int i = 0; i < q->dimensions[0]; i++) {
+      double q_value = *(double *)(q->data + i*q->strides[0]);
+      double phi_value = *(double *)(phi->data + i*phi->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      if (q_value == 0)
+          *result_value = 0.0;
+      else
+          *result_value = model->evaluate_rphi(q_value, phi_value);
+    }
+    return PyArray_Return(result); 
+ }
+ 
+ /**
+ * Function to call to evaluate model
+ * @param args: input numpy array  [x[],y[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDimXY( CylinderModel* model, 
+                              PyArrayObject *x, PyArrayObject *y)
+ {
+    PyArrayObject *result;
+    //check validity of input vectors
+    if (x->nd != 1 || x->descr->type_num != PyArray_DOUBLE
+        || y->nd != 1 || y->descr->type_num != PyArray_DOUBLE
+        || y->dimensions[0] != x->dimensions[0]){
+        
+        //const char * message= "Invalid array: x->nd=%d,type_num=%d\n",x->nd,x->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+	result= (PyArrayObject *)PyArray_FromDims(x->nd,(int*)(x->dimensions), PyArray_DOUBLE);
+
+	if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	}
+	
+    for (int i = 0; i < x->dimensions[0]; i++) {
+      double x_value = *(double *)(x->data + i*x->strides[0]);
+      double y_value = *(double *)(y->data + i*y->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      *result_value = (*model)(x_value, y_value);
+    }
+    return PyArray_Return(result); 
+ }
+ 
 /**
  * Function to call to evaluate model
@@ -194,21 +299,33 @@
 	    	return NULL;
 	    }
+	    
 	    // We have a vector q, get the q and phi values at which
 	    // to evaluate I(q,phi)
-	    q_value = CCylinderModel_readDouble(PyList_GET_ITEM(pars,0));
-	    phi_value = CCylinderModel_readDouble(PyList_GET_ITEM(pars,1));
-	    // Skip zero
-	    if (q_value==0) {
-	    	return Py_BuildValue("d",0.0);
-	    }
-		return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value));
-
+	    PyObject *q, *phi;
+	    q = PyList_GET_ITEM(pars,0);
+	    phi = PyList_GET_ITEM(pars,1);
+	    
+	    if (PyArray_Check(q) && PyArray_Check(phi)) {
+		  return evaluateTwoDim(self->model, (PyArrayObject*)q, (PyArrayObject*)phi); 
+	    } else if (PyArray_Check(q) || PyArray_Check(phi)) {
+	        return NULL;
+	    } else {
+	        q_value = CCylinderModel_readDouble(PyList_GET_ITEM(pars,0));
+	        phi_value = CCylinderModel_readDouble(PyList_GET_ITEM(pars,1));
+	        // Skip zero
+	        if (q_value==0) {
+	    	    return Py_BuildValue("d",0.0);
+	        }
+		    return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value));
+        }
 	} else {
-
-		// We have a scalar q, we will evaluate I(q)
+	if (PyArray_Check(pars)) {
+		return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+    } else {
+	    // We have a scalar q, we will evaluate I(q)
 		q_value = CCylinderModel_readDouble(pars);		
-		
-		return Py_BuildValue("d",(*(self->model))(q_value));
-	}	
+	    return Py_BuildValue("d",(*(self->model))(q_value));
+	    }	
+	}
 }
 
@@ -265,14 +382,28 @@
 	    // We have a vector q, get the qx and qy values at which
 	    // to evaluate I(qx,qy)
-	    qx_value = CCylinderModel_readDouble(PyList_GET_ITEM(pars,0));
-	    qy_value = CCylinderModel_readDouble(PyList_GET_ITEM(pars,1));
-	    return Py_BuildValue("d",(*(self->model))(qx_value,qy_value));
+	    PyObject *x, *y;
+	    x = PyList_GET_ITEM(pars,0);
+	    y = PyList_GET_ITEM(pars,1);
+	    
+	    if (PyArray_Check(x) && PyArray_Check(y)) {
+	       
+		    return evaluateTwoDimXY(self->model, (PyArrayObject*)x, (PyArrayObject*)y); 
+	    } else if (PyArray_Check(x) || PyArray_Check(y)) {
+	        return NULL;
+	    } else {    
+	        qx_value = CCylinderModel_readDouble(PyList_GET_ITEM(pars,0));
+	        qy_value = CCylinderModel_readDouble(PyList_GET_ITEM(pars,1));
+	        return Py_BuildValue("d",(*(self->model))(qx_value,qy_value));
+	    }
 
 	} else {
-
-		// We have a scalar q, we will evaluate I(q)
-		qx_value = CCylinderModel_readDouble(pars);		
+		if (PyArray_Check(pars)) {
+		    return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+        } else {
+		    // We have a scalar q, we will evaluate I(q)
+		    qx_value = CCylinderModel_readDouble(pars);		
 		
-		return Py_BuildValue("d",(*(self->model))(qx_value));
+		    return Py_BuildValue("d",(*(self->model))(qx_value));
+		}
 	}	
 }
@@ -376,8 +507,4 @@
 
 
-static PyMethodDef module_methods[] = {
-    {NULL} 
-};
-
 /**
  * Function used to add the model class to a module
Index: sansmodels/src/sans/models/c_models/CSphereModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CSphereModel.cpp	(revision 836fe6e73d15cf2f32c74bbf3b953b37d776f2ff)
+++ sansmodels/src/sans/models/c_models/CSphereModel.cpp	(revision 812b901e8b8451f555bc0649e5ee190827ff9130)
@@ -23,6 +23,10 @@
  */
  
+#define NO_IMPORT_ARRAY
+#define PY_ARRAY_UNIQUE_SYMBOL PyArray_API_sans
+
 extern "C" {
 #include <Python.h>
+#include <arrayobject.h>
 #include "structmember.h"
 #include <stdio.h>
@@ -130,6 +134,110 @@
     }
 }
-
-
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array q[] 
+ * @return: numpy array object 
+ */
+ 
+static PyObject *evaluateOneDim(SphereModel* model, PyArrayObject *q){
+    PyArrayObject *result;
+   
+    // Check validity of array q , q must be of dimension 1, an array of double
+    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE)
+    {
+        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message);
+        return NULL;
+    }
+    result = (PyArrayObject *)PyArray_FromDims(q->nd, (int *)(q->dimensions), 
+										  PyArray_DOUBLE);
+	if (result == NULL) {
+        const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+		return NULL;
+	}
+	 for (int i = 0; i < q->dimensions[0]; i++){
+      double q_value  = *(double *)(q->data + i*q->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      *result_value =(*model)(q_value);
+	}
+    return PyArray_Return(result); 
+ }
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array  [q[],phi[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDim( SphereModel* model, 
+                              PyArrayObject *q, PyArrayObject *phi)
+ {
+    PyArrayObject *result;
+    //check validity of input vectors
+    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE
+        || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE
+        || phi->dimensions[0] != q->dimensions[0]){
+        
+        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+	result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE);
+
+	if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	}
+	
+    for (int i = 0; i < q->dimensions[0]; i++) {
+      double q_value = *(double *)(q->data + i*q->strides[0]);
+      double phi_value = *(double *)(phi->data + i*phi->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      if (q_value == 0)
+          *result_value = 0.0;
+      else
+          *result_value = model->evaluate_rphi(q_value, phi_value);
+    }
+    return PyArray_Return(result); 
+ }
+ 
+ /**
+ * Function to call to evaluate model
+ * @param args: input numpy array  [x[],y[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDimXY( SphereModel* model, 
+                              PyArrayObject *x, PyArrayObject *y)
+ {
+    PyArrayObject *result;
+    //check validity of input vectors
+    if (x->nd != 1 || x->descr->type_num != PyArray_DOUBLE
+        || y->nd != 1 || y->descr->type_num != PyArray_DOUBLE
+        || y->dimensions[0] != x->dimensions[0]){
+        
+        //const char * message= "Invalid array: x->nd=%d,type_num=%d\n",x->nd,x->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+	result= (PyArrayObject *)PyArray_FromDims(x->nd,(int*)(x->dimensions), PyArray_DOUBLE);
+
+	if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	}
+	
+    for (int i = 0; i < x->dimensions[0]; i++) {
+      double x_value = *(double *)(x->data + i*x->strides[0]);
+      double y_value = *(double *)(y->data + i*y->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      if (x_value == 0)
+          *result_value = 0.0;
+      else
+          *result_value = (*model)(x_value, y_value);
+    }
+    return PyArray_Return(result); 
+ }
+ 
 /**
  * Function to call to evaluate model
@@ -173,21 +281,33 @@
 	    	return NULL;
 	    }
+	    
 	    // We have a vector q, get the q and phi values at which
 	    // to evaluate I(q,phi)
-	    q_value = CSphereModel_readDouble(PyList_GET_ITEM(pars,0));
-	    phi_value = CSphereModel_readDouble(PyList_GET_ITEM(pars,1));
-	    // Skip zero
-	    if (q_value==0) {
-	    	return Py_BuildValue("d",0.0);
-	    }
-		return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value));
-
+	    PyObject *q, *phi;
+	    q = PyList_GET_ITEM(pars,0);
+	    phi = PyList_GET_ITEM(pars,1);
+	    
+	    if (PyArray_Check(q) && PyArray_Check(phi)) {
+		  return evaluateTwoDim(self->model, (PyArrayObject*)q, (PyArrayObject*)phi); 
+	    } else if (PyArray_Check(q) || PyArray_Check(phi)) {
+	        return NULL;
+	    } else {
+	        q_value = CSphereModel_readDouble(PyList_GET_ITEM(pars,0));
+	        phi_value = CSphereModel_readDouble(PyList_GET_ITEM(pars,1));
+	        // Skip zero
+	        if (q_value==0) {
+	    	    return Py_BuildValue("d",0.0);
+	        }
+		    return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value));
+        }
 	} else {
-
-		// We have a scalar q, we will evaluate I(q)
+	if (PyArray_Check(pars)) {
+		return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+    } else {
+	    // We have a scalar q, we will evaluate I(q)
 		q_value = CSphereModel_readDouble(pars);		
-		
-		return Py_BuildValue("d",(*(self->model))(q_value));
-	}	
+	    return Py_BuildValue("d",(*(self->model))(q_value));
+	    }	
+	}
 }
 
@@ -235,14 +355,27 @@
 	    // We have a vector q, get the qx and qy values at which
 	    // to evaluate I(qx,qy)
-	    qx_value = CSphereModel_readDouble(PyList_GET_ITEM(pars,0));
-	    qy_value = CSphereModel_readDouble(PyList_GET_ITEM(pars,1));
-	    return Py_BuildValue("d",(*(self->model))(qx_value,qy_value));
+	    PyObject *x, *y;
+	    x = PyList_GET_ITEM(pars,0);
+	    y = PyList_GET_ITEM(pars,1);
+	    
+	    if (PyArray_Check(x) && PyArray_Check(y)) {
+		  return evaluateTwoDimXY(self->model, (PyArrayObject*)x, (PyArrayObject*)y); 
+	    } else if (PyArray_Check(x) || PyArray_Check(y)) {
+	        return NULL;
+	    } else {    
+	        qx_value = CSphereModel_readDouble(PyList_GET_ITEM(pars,0));
+	        qy_value = CSphereModel_readDouble(PyList_GET_ITEM(pars,1));
+	        return Py_BuildValue("d",(*(self->model))(qx_value,qy_value));
+	    }
 
 	} else {
-
-		// We have a scalar q, we will evaluate I(q)
-		qx_value = CSphereModel_readDouble(pars);		
+		if (PyArray_Check(pars)) {
+		    return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+        } else {
+		    // We have a scalar q, we will evaluate I(q)
+		    qx_value = CSphereModel_readDouble(pars);		
 		
-		return Py_BuildValue("d",(*(self->model))(qx_value));
+		    return Py_BuildValue("d",(*(self->model))(qx_value));
+		}
 	}	
 }
@@ -340,8 +473,4 @@
 
 
-static PyMethodDef module_methods[] = {
-    {NULL} 
-};
-
 /**
  * Function used to add the model class to a module
Index: sansmodels/src/sans/models/c_models/c_models.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/c_models.cpp	(revision 770e4c9e7d8f55e328b706d080c8771925265c1e)
+++ sansmodels/src/sans/models/c_models/c_models.cpp	(revision 812b901e8b8451f555bc0649e5ee190827ff9130)
@@ -7,4 +7,7 @@
 #include <Python.h>
 #include <parameters.hh>
+#define PY_ARRAY_UNIQUE_SYMBOL PyArray_API_sans
+#include "arrayobject.h"
+
 
 void addCCylinderModel(PyObject *module);
@@ -117,5 +120,5 @@
 	double *values;
 	double *weights;
-    int i;
+    //int i;
 
 	if (!PyArg_ParseTuple(args, "OOO", &disp, &val_obj, &wei_obj)) return NULL;
@@ -162,4 +165,5 @@
     m = Py_InitModule3("c_models", module_methods,
                        "C extension module for SANS scattering models.");
+	import_array();
 
 	addCCylinderModel(m);
Index: sansmodels/src/sans/models/c_models/classTemplate.txt
===================================================================
--- sansmodels/src/sans/models/c_models/classTemplate.txt	(revision af03ddd0932915d6659404ab7824e4223ea2944f)
+++ sansmodels/src/sans/models/c_models/classTemplate.txt	(revision 812b901e8b8451f555bc0649e5ee190827ff9130)
@@ -23,6 +23,10 @@
  */
  
+#define NO_IMPORT_ARRAY
+#define PY_ARRAY_UNIQUE_SYMBOL PyArray_API_sans
+
 extern "C" {
 #include <Python.h>
+#include <arrayobject.h>
 #include "structmember.h"
 #include <stdio.h>
@@ -116,6 +120,107 @@
     }
 }
-
-
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array q[] 
+ * @return: numpy array object 
+ */
+ 
+static PyObject *evaluateOneDim([CMODEL]* model, PyArrayObject *q){
+    PyArrayObject *result;
+   
+    // Check validity of array q , q must be of dimension 1, an array of double
+    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE)
+    {
+        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message);
+        return NULL;
+    }
+    result = (PyArrayObject *)PyArray_FromDims(q->nd, (int *)(q->dimensions), 
+										  PyArray_DOUBLE);
+	if (result == NULL) {
+        const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+		return NULL;
+	}
+	 for (int i = 0; i < q->dimensions[0]; i++){
+      double q_value  = *(double *)(q->data + i*q->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      *result_value =(*model)(q_value);
+	}
+    return PyArray_Return(result); 
+ }
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array  [q[],phi[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDim( [CMODEL]* model, 
+                              PyArrayObject *q, PyArrayObject *phi)
+ {
+    PyArrayObject *result;
+    //check validity of input vectors
+    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE
+        || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE
+        || phi->dimensions[0] != q->dimensions[0]){
+        
+        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+	result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE);
+
+	if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	}
+	
+    for (int i = 0; i < q->dimensions[0]; i++) {
+      double q_value = *(double *)(q->data + i*q->strides[0]);
+      double phi_value = *(double *)(phi->data + i*phi->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      if (q_value == 0)
+          *result_value = 0.0;
+      else
+          *result_value = model->evaluate_rphi(q_value, phi_value);
+    }
+    return PyArray_Return(result); 
+ }
+ 
+ /**
+ * Function to call to evaluate model
+ * @param args: input numpy array  [x[],y[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDimXY( [CMODEL]* model, 
+                              PyArrayObject *x, PyArrayObject *y)
+ {
+    PyArrayObject *result;
+    //check validity of input vectors
+    if (x->nd != 1 || x->descr->type_num != PyArray_DOUBLE
+        || y->nd != 1 || y->descr->type_num != PyArray_DOUBLE
+        || y->dimensions[0] != x->dimensions[0]){
+        
+        //const char * message= "Invalid array: x->nd=%d,type_num=%d\n",x->nd,x->descr->type_num;
+        //PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+	result= (PyArrayObject *)PyArray_FromDims(x->nd,(int*)(x->dimensions), PyArray_DOUBLE);
+
+	if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	}
+	
+    for (int i = 0; i < x->dimensions[0]; i++) {
+      double x_value = *(double *)(x->data + i*x->strides[0]);
+      double y_value = *(double *)(y->data + i*y->strides[0]);
+      double *result_value = (double *)(result->data + i*result->strides[0]);
+      *result_value = (*model)(x_value, y_value);
+    }
+    return PyArray_Return(result); 
+ }
+ 
 /**
  * Function to call to evaluate model
@@ -149,21 +254,33 @@
 	    	return NULL;
 	    }
+	    
 	    // We have a vector q, get the q and phi values at which
 	    // to evaluate I(q,phi)
-	    q_value = [PYTHONCLASS]_readDouble(PyList_GET_ITEM(pars,0));
-	    phi_value = [PYTHONCLASS]_readDouble(PyList_GET_ITEM(pars,1));
-	    // Skip zero
-	    if (q_value==0) {
-	    	return Py_BuildValue("d",0.0);
-	    }
-		return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value));
-
+	    PyObject *q, *phi;
+	    q = PyList_GET_ITEM(pars,0);
+	    phi = PyList_GET_ITEM(pars,1);
+	    
+	    if (PyArray_Check(q) && PyArray_Check(phi)) {
+		  return evaluateTwoDim(self->model, (PyArrayObject*)q, (PyArrayObject*)phi); 
+	    } else if (PyArray_Check(q) || PyArray_Check(phi)) {
+	        return NULL;
+	    } else {
+	        q_value = [PYTHONCLASS]_readDouble(PyList_GET_ITEM(pars,0));
+	        phi_value = [PYTHONCLASS]_readDouble(PyList_GET_ITEM(pars,1));
+	        // Skip zero
+	        if (q_value==0) {
+	    	    return Py_BuildValue("d",0.0);
+	        }
+		    return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value));
+        }
 	} else {
-
-		// We have a scalar q, we will evaluate I(q)
+	if (PyArray_Check(pars)) {
+		return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+    } else {
+	    // We have a scalar q, we will evaluate I(q)
 		q_value = [PYTHONCLASS]_readDouble(pars);		
-		
-		return Py_BuildValue("d",(*(self->model))(q_value));
-	}	
+	    return Py_BuildValue("d",(*(self->model))(q_value));
+	    }	
+	}
 }
 
@@ -201,14 +318,28 @@
 	    // We have a vector q, get the qx and qy values at which
 	    // to evaluate I(qx,qy)
-	    qx_value = [PYTHONCLASS]_readDouble(PyList_GET_ITEM(pars,0));
-	    qy_value = [PYTHONCLASS]_readDouble(PyList_GET_ITEM(pars,1));
-	    return Py_BuildValue("d",(*(self->model))(qx_value,qy_value));
+	    PyObject *x, *y;
+	    x = PyList_GET_ITEM(pars,0);
+	    y = PyList_GET_ITEM(pars,1);
+	    
+	    if (PyArray_Check(x) && PyArray_Check(y)) {
+	       
+		    return evaluateTwoDimXY(self->model, (PyArrayObject*)x, (PyArrayObject*)y); 
+	    } else if (PyArray_Check(x) || PyArray_Check(y)) {
+	        return NULL;
+	    } else {    
+	        qx_value = [PYTHONCLASS]_readDouble(PyList_GET_ITEM(pars,0));
+	        qy_value = [PYTHONCLASS]_readDouble(PyList_GET_ITEM(pars,1));
+	        return Py_BuildValue("d",(*(self->model))(qx_value,qy_value));
+	    }
 
 	} else {
-
-		// We have a scalar q, we will evaluate I(q)
-		qx_value = [PYTHONCLASS]_readDouble(pars);		
+		if (PyArray_Check(pars)) {
+		    return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+        } else {
+		    // We have a scalar q, we will evaluate I(q)
+		    qx_value = [PYTHONCLASS]_readDouble(pars);		
 		
-		return Py_BuildValue("d",(*(self->model))(qx_value));
+		    return Py_BuildValue("d",(*(self->model))(qx_value));
+		}
 	}	
 }
@@ -302,8 +433,4 @@
 
 
-static PyMethodDef module_methods[] = {
-    {NULL} 
-};
-
 /**
  * Function used to add the model class to a module
Index: sansmodels/src/sans/models/test/utest_models_array.py
===================================================================
--- sansmodels/src/sans/models/test/utest_models_array.py	(revision 812b901e8b8451f555bc0649e5ee190827ff9130)
+++ sansmodels/src/sans/models/test/utest_models_array.py	(revision 812b901e8b8451f555bc0649e5ee190827ff9130)
@@ -0,0 +1,94 @@
+"""
+    Unit tests for specific models using numpy array input
+    @author: Gervaise B Alina/ UTK
+"""
+
+import unittest, time, math, numpy
+
+        
+class TestSphere(unittest.TestCase):
+    """ Unit tests for sphere model """
+    
+    def setUp(self):
+        from sans.models.SphereModel import SphereModel
+        self.comp = SphereModel()
+        self.x = numpy.array([1.0,2.0,3.0, 4.0])
+        self.y = self.x +1
+        
+    def test1D(self):
+        """ Test 1D model for a sphere  with vector as input"""
+        answer=numpy.array([5.63877831e-05,2.57231782e-06,2.73704050e-07,2.54229069e-08])
+       
+        testvector= self.comp.run(self.x)
+       
+        self.assertAlmostEqual(len(testvector),4)
+        for i in xrange(len(answer)):
+            self.assertAlmostEqual(testvector[i],answer[i])
+       
+    def test1D_1(self):
+        """ Test 2D model for a sphere  with scalar as input"""
+        self.assertAlmostEqual(self.comp.run(1.0),5.63877831e-05, 4)
+         
+    def test1D_2(self):
+        """ Test 2D model for a sphere for 2 scalar """
+        self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 56.3878e-06, 4)
+   
+        
+class TestCylinder(unittest.TestCase):
+    """ Unit tests for sphere model """
+    
+    def setUp(self):
+        from sans.models.CylinderModel import CylinderModel
+        self.comp = CylinderModel()
+        self.x = numpy.array([1.0,2.0,3.0, 4.0])
+        self.y = numpy.array([1.0,2.0,3.0, 4.0])
+        
+    def test1D(self):
+        """ Test 1D model for a Cylinder  with vector as input"""
+        answer = numpy.array([1.98860592e-04,7.03686335e-05,2.89144683e-06,2.04282827e-06])
+        testvector= self.comp.run(self.x)
+       
+        self.assertAlmostEqual(len(testvector),4)
+        for i in xrange(len(answer)):
+            self.assertAlmostEqual(testvector[i],answer[i])
+       
+    def test1D_1(self):
+        """ Test 2D model for a Cylinder with scalar as input"""
+        self.assertAlmostEqual(self.comp.run(1.0),1.9886e-04, 4)
+         
+    def test1D_2(self):
+        """ Test 2D model for a Cylinder for 2 scalar """
+        self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 56.3878e-06, 4)
+        
+    def test1D_3(self):
+        """ Test 2D model for a Cylinder for 2 vectors as input """
+        ans_input = numpy.zeros(len(self.x))
+        temp_x = numpy.zeros(len(self.x))
+        temp_y = numpy.zeros(len(self.y))
+        
+        for i in xrange(len(self.x)):
+            qx = self.x[i]
+            qy = self.y[i]
+          
+            temp_x[i]= qx*math.cos(qy)
+            temp_y[i]= qx*math.sin(qy)
+            
+            value = math.sqrt(temp_x[i]*temp_x[i]+ temp_y[i]*temp_y[i] )#qx*qx +qy*qy)
+            ans_input[i]= value
+         
+        vect_runXY_qx_qy = self.comp.runXY([temp_x, temp_y])
+        vect_run_x_y = self.comp.run([self.x, self.y])
+        
+        for i in xrange(len(vect_runXY_qx_qy)):
+            self.assertAlmostEqual(vect_runXY_qx_qy[i], vect_run_x_y[i])
+            
+        vect_run_x = self.comp.run(self.x)
+        vect_run_answer = self.comp.run(ans_input)
+        
+        for i in xrange(len(vect_run_x )):
+            self.assertAlmostEqual(vect_run_x [i], vect_run_answer[i])
+            
+
+
+if __name__ == '__main__':
+    unittest.main()
