Index: sansmodels/src/sans/models/c_models/CCylinder.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CCylinder.cpp	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/CCylinder.cpp	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,354 @@
+/** CCylinder
+ *
+ * C extension
+ *
+ */
+extern "C" {
+#include <Python.h>
+#include "structmember.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "cylinder.h"
+}
+
+
+#include "models.hh"
+#include "dispersion_visitor.hh"
+
+/// Error object for raised exceptions
+static PyObject * CCylinderModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    PyObject * dispersion;
+    Cylinder * model;
+    /// Log for unit testing
+    PyObject * log;
+    /// Model parameters
+	CylinderParameters model_pars;
+} CCylinderModel;
+
+
+static void
+CCylinderModel_dealloc(CCylinderModel* self)
+{
+    self->ob_type->tp_free((PyObject*)self);
+
+
+}
+
+static PyObject *
+CCylinderModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CCylinderModel *self;
+
+    self = (CCylinderModel *)type->tp_alloc(type, 0);
+
+    return (PyObject *)self;
+}
+
+static int
+CCylinderModel_init(CCylinderModel *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	printf("Hello cylinder\n");
+    	// Create parameters
+        self->params = PyDict_New();
+        self->dispersion = PyDict_New();
+        self->model = new Cylinder();
+
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"scale",     Py_BuildValue("d", self->model->scale()));
+        PyDict_SetItemString(self->params,"length",    Py_BuildValue("d", self->model->length()));
+        PyDict_SetItemString(self->params,"cyl_theta", Py_BuildValue("d", self->model->cyl_theta()));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d", self->model->background()));
+        PyDict_SetItemString(self->params,"radius",    Py_BuildValue("d", self->model->radius()));
+        PyDict_SetItemString(self->params,"contrast",  Py_BuildValue("d", self->model->contrast()));
+        PyDict_SetItemString(self->params,"cyl_phi",   Py_BuildValue("d", self->model->cyl_phi()));
+
+        // Initialize dispersion / averaging parameter dict
+        DispersionVisitor* visitor = new DispersionVisitor();
+        PyObject * disp_dict = PyDict_New();
+        self->model->radius.dispersion->accept_as_source(visitor, self->model->radius.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "disp_radius", disp_dict);
+
+        disp_dict = PyDict_New();
+        self->model->length.dispersion->accept_as_source(visitor, self->model->length.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "disp_length", disp_dict);
+
+        disp_dict = PyDict_New();
+        self->model->cyl_phi.dispersion->accept_as_source(visitor, self->model->cyl_phi.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "disp_cyl_phi", disp_dict);
+
+        disp_dict = PyDict_New();
+        self->model->cyl_theta.dispersion->accept_as_source(visitor, self->model->cyl_theta.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "disp_cyl_theta", disp_dict);
+
+        // Create empty log
+        self->log = PyDict_New();
+
+    }
+    return 0;
+}
+
+static PyMemberDef CCylinderModel_members[] = {
+	{"params", T_OBJECT, offsetof(CCylinderModel, params), 0,
+	 "Parameters"},
+	{"dispersion", T_OBJECT, offsetof(CCylinderModel, dispersion), 0,
+	  "Dispersion parameters"},
+    {"log", T_OBJECT, offsetof(CCylinderModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CCylinderModel_readDouble(PyObject *p) {
+    if (PyFloat_Check(p)==1) {
+        return (double)(((PyFloatObject *)(p))->ob_fval);
+    } else if (PyInt_Check(p)==1) {
+        return (double)(((PyIntObject *)(p))->ob_ival);
+    } else if (PyLong_Check(p)==1) {
+        return (double)PyLong_AsLong(p);
+    } else {
+        return 0.0;
+    }
+}
+
+/**
+ * Method to set the dispersion model for a parameter.
+ * We need to update the dispersion object of the parameter
+ * on the C++ side and update the dispersion dictionary on the python side.
+ *
+ * TODO: define read-only integers that will be data members of this
+ * class and allow the user to specify which model is needed with a code name.
+ *
+ * This method should take in a code name (FLAT, GAUSSIAN) and a string
+ * that represent the parameter to attach the dispersion model to.
+ *
+ */
+static PyObject * set_dispersion_model(CCylinderModel *self, PyObject *args) {
+
+}
+
+
+/**
+ * Function to call to evaluate model
+ * @param args: input q or [q,phi]
+ * @return: function value
+ */
+static PyObject * run(CCylinderModel *self, PyObject *args) {
+	double q_value, phi_value;
+	PyObject* pars;
+	int npars;
+
+	// Get parameters
+
+	// Reader parameter dictionary
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model->cyl_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_theta") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") );
+    self->model->cyl_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_phi") );
+
+    // Read in dispersion parameters
+    DispersionVisitor* visitor = new DispersionVisitor();
+    PyObject* disp_dict = PyDict_New();
+    disp_dict = PyDict_GetItemString(self->dispersion, "disp_radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict);
+    self->model->cyl_theta.dispersion->accept_as_destination(visitor, self->model->cyl_theta.dispersion, disp_dict);
+    self->model->cyl_phi.dispersion->accept_as_destination(visitor, self->model->cyl_phi.dispersion, disp_dict);
+
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CCylinderModelError,
+	    	"CCylinderModel.run expects a q value.");
+		return NULL;
+	}
+
+	// Check params
+	if( PyList_Check(pars)==1) {
+
+		// Length of list should be 2 for I(q,phi)
+	    npars = PyList_GET_SIZE(pars);
+	    if(npars!=2) {
+	    	PyErr_SetString(CCylinderModelError,
+	    		"CCylinderModel.run expects a double or a list of dimension 2.");
+	    	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));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CCylinderModel_readDouble(pars);
+
+		return Py_BuildValue("d",(*(self->model))(q_value));
+	}
+}
+
+/**
+ * Function to call to evaluate model in cartesian coordinates
+ * @param args: input q or [qx, qy]]
+ * @return: function value
+ */
+static PyObject * runXY(CCylinderModel *self, PyObject *args) {
+	double qx_value, qy_value;
+	PyObject* pars;
+	int npars;
+
+	// Get parameters
+
+	// Reader parameter dictionary
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model->cyl_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_theta") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") );
+    self->model->cyl_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_phi") );
+
+    // Read in dispersion parameters
+    DispersionVisitor* visitor = new DispersionVisitor();
+    PyObject* disp_dict = PyDict_New();
+    disp_dict = PyDict_GetItemString(self->dispersion, "disp_radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict);
+    self->model->cyl_theta.dispersion->accept_as_destination(visitor, self->model->cyl_theta.dispersion, disp_dict);
+    self->model->cyl_phi.dispersion->accept_as_destination(visitor, self->model->cyl_phi.dispersion, disp_dict);
+
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CCylinderModelError,
+	    	"CCylinderModel.run expects a q value.");
+		return NULL;
+	}
+
+	// Check params
+	if( PyList_Check(pars)==1) {
+
+		// Length of list should be 2 for I(qx, qy))
+	    npars = PyList_GET_SIZE(pars);
+	    if(npars!=2) {
+	    	PyErr_SetString(CCylinderModelError,
+	    		"CCylinderModel.run expects a double or a list of dimension 2.");
+	    	return NULL;
+	    }
+	    // 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));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		qx_value = CCylinderModel_readDouble(pars);
+
+		return Py_BuildValue("d",(*(self->model))(qx_value));
+	}
+}
+
+static PyObject * reset(CCylinderModel *self, PyObject *args) {
+
+
+    return Py_BuildValue("d",0.0);
+}
+
+
+static PyMethodDef CCylinderModel_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+    {"runXY",      (PyCFunction)runXY     , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy"},
+    {"reset",    (PyCFunction)reset   , METH_VARARGS,
+      "Reset pair correlation"},
+    //{"numerical_1D",      (PyCFunction)numerical_1D     , METH_VARARGS,
+    //  "Evaluate the 1D model at a given Q"},
+   {NULL}
+};
+
+static PyTypeObject CCylinderModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CCylinderModel",             /*tp_name*/
+    sizeof(CCylinderModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CCylinderModel_dealloc, /*tp_dealloc*/
+    0,                         /*tp_print*/
+    0,                         /*tp_getattr*/
+    0,                         /*tp_setattr*/
+    0,                         /*tp_compare*/
+    0,                         /*tp_repr*/
+    0,                         /*tp_as_number*/
+    0,                         /*tp_as_sequence*/
+    0,                         /*tp_as_mapping*/
+    0,                         /*tp_hash */
+    0,                         /*tp_call*/
+    0,                         /*tp_str*/
+    0,                         /*tp_getattro*/
+    0,                         /*tp_setattro*/
+    0,                         /*tp_as_buffer*/
+    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
+    "CCylinderModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CCylinderModel_methods,             /* tp_methods */
+    CCylinderModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CCylinderModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CCylinderModel_new,                 /* tp_new */
+};
+
+
+static PyMethodDef module_methods[] = {
+    {NULL}
+};
+
+/**
+ * Function used to add the model class to a module
+ * @param module: module to add the class to
+ */
+void addCCylinderModel(PyObject *module) {
+	PyObject *d;
+
+    if (PyType_Ready(&CCylinderModelType) < 0)
+        return;
+
+    Py_INCREF(&CCylinderModelType);
+    PyModule_AddObject(module, "CCylinderModel", (PyObject *)&CCylinderModelType);
+
+    d = PyModule_GetDict(module);
+    CCylinderModelError = PyErr_NewException("CCylinderModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CCylinderModelError", CCylinderModelError);
+}
+
Index: sansmodels/src/sans/models/c_models/cylinder.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/cylinder.cpp	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/cylinder.cpp	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,233 @@
+/**
+	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 2008, University of Tennessee
+ */
+
+/**
+ * Scattering model classes
+ * The classes use the IGOR library found in
+ *   sansmodels/src/libigor
+ *
+ *	TODO: refactor so that we pull in the old sansmodels.c_extensions
+ */
+
+#include <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "libCylinder.h"
+	#include "cylinder.h"
+}
+
+Cylinder :: Cylinder() {
+	scale      = Parameter(1.0);
+	radius     = Parameter(20.0, true);
+	radius.set_min(0.0);
+	length     = Parameter(400.0, true);
+	length.set_min(0.0);
+	contrast   = Parameter(3.e-6);
+	background = Parameter(0.0);
+	cyl_theta  = Parameter(0.0, true);
+	cyl_phi    = Parameter(0.0, true);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * The NIST IGOR library is used for the actual calculation.
+ * @param q: q-value
+ * @return: function value
+ */
+double Cylinder :: operator()(double q) {
+	double dp[5];
+
+	// Fill parameter array for IGOR library
+	// Add the background after averaging
+	dp[0] = scale();
+	dp[1] = radius();
+	dp[2] = length();
+	dp[3] = contrast();
+	dp[4] = 0.0;
+
+	// Get the dispersion points for the radius
+	vector<WeightPoint> weights_rad;
+	radius.get_weights(weights_rad);
+
+	// Get the dispersion points for the length
+	vector<WeightPoint> weights_len;
+	length.get_weights(weights_len);
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+
+	// Loop over radius weight points
+	for(int i=0; i<weights_rad.size(); i++) {
+		dp[1] = weights_rad[i].value;
+
+		// Loop over length weight points
+		for(int j=0; j<weights_len.size(); j++) {
+			dp[2] = weights_len[j].value;
+
+			sum += weights_rad[i].weight
+				* weights_len[j].weight * CylinderForm(dp, q);
+			norm += weights_rad[i].weight
+				* weights_len[j].weight;
+		}
+	}
+	return sum/norm + background();
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double Cylinder :: operator()(double qx, double qy) {
+	CylinderParameters dp;
+	// Fill parameter array
+	dp.scale      = scale();
+	dp.radius     = radius();
+	dp.length     = length();
+	dp.contrast   = contrast();
+	dp.background = 0.0;
+	dp.cyl_theta  = cyl_theta();
+	dp.cyl_phi    = cyl_phi();
+
+	// Get the dispersion points for the radius
+	vector<WeightPoint> weights_rad;
+	radius.get_weights(weights_rad);
+
+	// Get the dispersion points for the length
+	vector<WeightPoint> weights_len;
+	length.get_weights(weights_len);
+
+	// Get angular averaging for theta
+	vector<WeightPoint> weights_theta;
+	cyl_theta.get_weights(weights_theta);
+
+	// Get angular averaging for phi
+	vector<WeightPoint> weights_phi;
+	cyl_phi.get_weights(weights_phi);
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+
+	// Loop over radius weight points
+	for(int i=0; i<weights_rad.size(); i++) {
+		dp.radius = weights_rad[i].value;
+
+
+		// Loop over length weight points
+		for(int j=0; j<weights_len.size(); j++) {
+			dp.length = weights_len[j].value;
+
+			// Average over theta distribution
+			for(int k=0; k<weights_theta.size(); k++) {
+				dp.cyl_theta = weights_theta[k].value;
+
+				// Average over phi distribution
+				for(int l=0; l<weights_phi.size(); l++) {
+					dp.cyl_phi = weights_phi[l].value;
+
+					double _ptvalue = weights_rad[i].weight
+						* weights_len[j].weight
+						* weights_theta[k].weight
+						* weights_phi[l].weight
+						* cylinder_analytical_2DXY(&dp, qx, qy);
+					if (weights_theta.size()>1) {
+						_ptvalue *= sin(weights_theta[k].value);
+					}
+					sum += _ptvalue;
+
+					norm += weights_rad[i].weight
+						* weights_len[j].weight
+						* weights_theta[k].weight
+						* weights_phi[l].weight;
+
+				}
+			}
+		}
+	}
+	// Averaging in theta needs an extra normalization
+	// factor to account for the sin(theta) term in the
+	// integration (see documentation).
+	if (weights_theta.size()>1) norm = norm / asin(1.0);
+	return sum/norm + background();
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double Cylinder :: evaluate_rphi(double q, double phi) {
+	double qx = q*cos(phi);
+	double qy = q*sin(phi);
+	return (*this).operator()(qx, qy);
+}
+
+// Testing code
+int main(void)
+{
+	Cylinder c = Cylinder();
+
+	printf("Length = %g\n", c.length());
+	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
+	c.radius.dispersion = new GaussianDispersion();
+	c.radius.dispersion->npts = 100;
+	c.radius.dispersion->width = 5;
+
+	//c.length.dispersion = GaussianDispersion();
+	//c.length.dispersion.npts = 20;
+	//c.length.dispersion.width = 65;
+
+	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
+	c.scale = 10.0;
+	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
+	printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
+	printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854));
+
+	// Average over phi at theta=90 deg
+	c.cyl_theta = 1.57;
+	double values_th[100];
+	double values[100];
+	double weights[100];
+	double pi = acos(-1.0);
+	printf("pi=%g\n", pi);
+	for(int i=0; i<100; i++){
+		values[i] = (float)i*2.0*pi/99.0;
+		values_th[i] = (float)i*pi/99.0;
+		weights[i] = 1.0;
+	}
+	//c.radius.dispersion->width = 0;
+	c.cyl_phi.dispersion = new ArrayDispersion();
+	c.cyl_theta.dispersion = new ArrayDispersion();
+	(*c.cyl_phi.dispersion).set_weights(100, values, weights);
+	(*c.cyl_theta.dispersion).set_weights(100, values_th, weights);
+
+	double i_avg = c(0.01, 0.01);
+	double i_1d = c(sqrt(0.0002));
+
+	printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg);
+	printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d);
+	printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg);
+
+
+	return 0;
+}
Index: sansmodels/src/sans/models/c_models/dispersion_visitor.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/dispersion_visitor.cpp	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/dispersion_visitor.cpp	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,63 @@
+#include "dispersion_visitor.hh"
+#include "parameters.hh"
+
+#ifndef __MODELS_STANDALONE__
+extern "C" {
+	#include <Python.h>
+}
+#endif
+
+void DispersionVisitor:: dispersion_to_dict(void* dispersion, void* dictionary) {
+#ifndef __MODELS_STANDALONE__
+	DispersionModel * disp = (DispersionModel*)dispersion;
+	PyObject * dict = (PyObject*)dictionary;
+
+	PyDict_SetItemString(dict, "type",  Py_BuildValue("s", "flat"));
+    PyDict_SetItemString(dict, "npts",  Py_BuildValue("i", disp->npts));
+    PyDict_SetItemString(dict, "width", Py_BuildValue("d", disp->width));
+#endif
+}
+
+void DispersionVisitor:: gaussian_to_dict(void* dispersion, void* dictionary) {
+#ifndef __MODELS_STANDALONE__
+	GaussianDispersion * disp = (GaussianDispersion*)dispersion;
+	PyObject * dict = (PyObject*)dictionary;
+
+	PyDict_SetItemString(dict, "type",  Py_BuildValue("s", "gaussian"));
+    PyDict_SetItemString(dict, "npts",  Py_BuildValue("i", disp->npts));
+    PyDict_SetItemString(dict, "width", Py_BuildValue("d", disp->width));
+    PyDict_SetItemString(dict, "nsigmas", Py_BuildValue("i", disp->nsigmas));
+#endif
+}
+
+void DispersionVisitor:: array_to_dict(void* dispersion, void* dictionary) {
+#ifndef __MODELS_STANDALONE__
+	ArrayDispersion * disp = (ArrayDispersion*)dispersion;
+	PyObject * dict = (PyObject*)dictionary;
+
+    PyDict_SetItemString(dict, "type",  Py_BuildValue("s", "array"));
+#endif
+}
+
+void DispersionVisitor:: dispersion_from_dict(void* dispersion, void* dictionary) {
+#ifndef __MODELS_STANDALONE__
+	DispersionModel * disp = (DispersionModel*)dispersion;
+	PyObject * dict = (PyObject*)dictionary;
+
+	disp->npts  = PyInt_AsLong( PyDict_GetItemString(dict, "npts") );
+	disp->width = PyFloat_AsDouble( PyDict_GetItemString(dict, "width") );
+#endif
+}
+
+void DispersionVisitor:: gaussian_from_dict(void* dispersion, void* dictionary) {
+#ifndef __MODELS_STANDALONE__
+	GaussianDispersion * disp = (GaussianDispersion*)dispersion;
+	PyObject * dict = (PyObject*)dictionary;
+
+	disp->npts    = PyInt_AsLong( PyDict_GetItemString(dict, "npts") );
+	disp->width   = PyFloat_AsDouble( PyDict_GetItemString(dict, "width") );
+	disp->nsigmas = PyFloat_AsDouble( PyDict_GetItemString(dict, "nsigmas") );
+#endif
+}
+
+void DispersionVisitor:: array_from_dict(void* dispersion, void* dictionary) {}
Index: sansmodels/src/sans/models/c_models/dispersion_visitor.hh
===================================================================
--- sansmodels/src/sans/models/c_models/dispersion_visitor.hh	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/dispersion_visitor.hh	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,30 @@
+/**
+	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 2008, University of Tennessee
+ */
+#ifndef DISP_VISIT_H
+#define DISP_VISIT_H
+
+class DispersionVisitor {
+public:
+	void dispersion_to_dict(void *, void *);
+	void gaussian_to_dict(void *, void *);
+	void array_to_dict(void *, void *);
+
+	void dispersion_from_dict(void*, void *);
+	void gaussian_from_dict(void*, void *);
+	void array_from_dict(void*, void *);
+
+};
+
+
+#endif
Index: sansmodels/src/sans/models/c_models/models.hh
===================================================================
--- sansmodels/src/sans/models/c_models/models.hh	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/models.hh	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,47 @@
+/**
+	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 2008, University of Tennessee
+ */
+#ifndef MODEL_CLASS_H
+#define MODEL_CLASS_H
+
+#include <vector>
+#include "parameters.hh"
+extern "C" {
+	#include "cylinder.h"
+}
+
+using namespace std;
+
+class Cylinder{
+
+public:
+	// Model parameters
+	Parameter radius;
+	Parameter scale;
+	Parameter length;
+	Parameter contrast;
+	Parameter background;
+	Parameter cyl_theta;
+	Parameter cyl_phi;
+	// TODO: replace this by an array of parameters
+
+	// Constructor
+	Cylinder();
+
+	// Operators to get I(Q)
+	double operator()(double q);
+	double operator()(double qx, double qy);
+	double evaluate_rphi(double q, double phi);
+};
+
+#endif
Index: sansmodels/src/sans/models/c_models/notes.txt
===================================================================
--- sansmodels/src/sans/models/c_models/notes.txt	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/notes.txt	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,13 @@
+The design of the SANS models is not easily scalable. Although a more efficient
+way (in terms of design and speed) as been implemented on the C/C++ side, the python
+design should be revisited. Areas of improvements include:
+
+1- [C side] Incorporate the C functions that use libigor directly into the new C++ class.
+  	That would mean rewritting the script generator to eliminate all the .c files.
+  
+2- Re-design the sans.models.BaseComponent class to store the parameters as objects,
+	like it's done on the C++ side. A better design may allow us to avoid the awful
+	hard-coding of the translation between the python dictionary holding the 
+	parameters and the C++ Parameter class.
+	 
+ 
Index: sansmodels/src/sans/models/c_models/parameters.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/parameters.cpp	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/parameters.cpp	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,231 @@
+/**
+	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 2008, University of Tennessee
+ */
+#include "parameters.hh"
+#include <stdio.h>
+#include <math.h>
+using namespace std;
+
+/**
+ * TODO: normalize all dispersion weight lists
+ */
+
+
+/**
+ * Weight points
+ */
+WeightPoint :: WeightPoint() {
+	value = 0.0;
+	weight = 0.0;
+}
+
+WeightPoint :: WeightPoint(double v, double w) {
+	value = v;
+	weight = w;
+}
+
+/**
+ * Dispersion models
+ */
+DispersionModel :: DispersionModel() {
+	npts  = 1;
+	width = 0.0;
+};
+
+void DispersionModel :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
+	visitor->dispersion_to_dict(from, to);
+}
+void DispersionModel :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
+	visitor->dispersion_from_dict(from, to);
+}
+void DispersionModel :: operator() (void *param, vector<WeightPoint> &weights){
+	// Check against zero width
+	if (width<=0) {
+		width = 0.0;
+		npts  = 1;
+	}
+
+	Parameter* par = (Parameter*)param;
+	double value = (*par)();
+
+	if (npts<2) {
+		weights.insert(weights.end(), WeightPoint(value, 1.0));
+	} else {
+		for(int i=0; i<npts; i++) {
+			double val = value + width * (1.0*i/float(npts-1) - 0.5);
+
+			if ( ((*par).has_min==false || val>(*par).min)
+			  && ((*par).has_max==false || val<(*par).max)  )
+				weights.insert(weights.end(), WeightPoint(val, 1.0));
+		}
+	}
+}
+
+/**
+ * Method to set the weights
+ * Not implemented for this class
+ */
+void DispersionModel :: set_weights(int npoints, double* values, double* weights){}
+
+/**
+ * Gaussian dispersion
+ */
+
+GaussianDispersion :: GaussianDispersion() {
+	npts  = 1;
+	width = 0.0;
+	nsigmas = 2;
+};
+
+void GaussianDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
+	visitor->gaussian_to_dict(from, to);
+}
+void GaussianDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
+	visitor->gaussian_from_dict(from, to);
+}
+
+double gaussian_weight(double mean, double sigma, double x) {
+	double vary, expo_value;
+    vary = x-mean;
+    expo_value = -vary*vary/(2*sigma*sigma);
+    //return 1.0;
+    return exp(expo_value);
+}
+
+/**
+ * Gaussian dispersion
+ * @param mean: mean value of the Gaussian
+ * @param sigma: standard deviation of the Gaussian
+ * @param x: value at which the Gaussian is evaluated
+ * @return: value of the Gaussian
+ */
+void GaussianDispersion :: operator() (void *param, vector<WeightPoint> &weights){
+	// Check against zero width
+	if (width<=0) {
+		width = 0.0;
+		npts  = 1;
+	}
+
+	Parameter* par = (Parameter*)param;
+	double value = (*par)();
+
+	if (npts<2) {
+		weights.insert(weights.end(), WeightPoint(value, 1.0));
+	} else {
+		for(int i=0; i<npts; i++) {
+			// We cover 2 sigmas on each side of the mean
+			double val = value + width * (4.0*i/float(npts-1) - 2.0);
+
+			if ( ((*par).has_min==false || val>(*par).min)
+			  && ((*par).has_max==false || val<(*par).max)  ) {
+				double _w = gaussian_weight(value, width, val);
+				weights.insert(weights.end(), WeightPoint(val, _w));
+			}
+		}
+	}
+}
+
+/**
+ * Array dispersion based on input arrays
+ */
+
+void ArrayDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
+	visitor->array_to_dict(from, to);
+}
+void ArrayDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
+	visitor->array_from_dict(from, to);
+}
+
+/**
+ * Method to get the weights
+ */
+void ArrayDispersion :: operator() (void *param, vector<WeightPoint> &weights) {
+	Parameter* par = (Parameter*)param;
+	double value = (*par)();
+
+	if (npts<2) {
+		weights.insert(weights.end(), WeightPoint(value, 1.0));
+	} else {
+		for(int i=0; i<npts; i++) {
+			if ( ((*par).has_min==false || _values[i]>(*par).min)
+			  && ((*par).has_max==false || _values[i]<(*par).max)  )
+				weights.insert(weights.end(), WeightPoint(_values[i], _weights[i]));
+		}
+	}
+}
+
+
+/**
+ * Method to set the weights
+ */
+void ArrayDispersion :: set_weights(int npoints, double* values, double* weights){
+	npts = npoints;
+	_values = values;
+	_weights = weights;
+}
+
+
+/**
+ * Parameters
+ */
+Parameter :: Parameter() {
+	value = 0;
+	min   = 0.0;
+	max   = 0.0;
+	has_min = false;
+	has_max = false;
+	has_dispersion = false;
+	dispersion = new GaussianDispersion();
+}
+
+Parameter :: Parameter(double _value) {
+	value = _value;
+	min   = 0.0;
+	max   = 0.0;
+	has_min = false;
+	has_max = false;
+	has_dispersion = false;
+	dispersion = new GaussianDispersion();
+}
+
+Parameter :: Parameter(double _value, bool disp) {
+	value = _value;
+	min   = 0.0;
+	max   = 0.0;
+	has_min = false;
+	has_max = false;
+	has_dispersion = disp;
+	dispersion = new GaussianDispersion();
+}
+
+void Parameter :: get_weights(vector<WeightPoint> &weights) {
+	(*dispersion)((void*)this, weights);
+}
+
+void Parameter :: set_min(double value) {
+	has_min = true;
+	min = value;
+}
+
+void Parameter :: set_max(double value) {
+	has_max = true;
+	max = value;
+}
+
+double Parameter :: operator()() {
+	return value;
+}
+
+double Parameter :: operator=(double _value){
+	value = _value;
+}
Index: sansmodels/src/sans/models/c_models/parameters.hh
===================================================================
--- sansmodels/src/sans/models/c_models/parameters.hh	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/parameters.hh	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,133 @@
+/**
+	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 2008, University of Tennessee
+ */
+#ifndef PARAM_CLASS_H
+#define PARAM_CLASS_H
+/**
+ * TODO: will need to write a bridge class
+ * to convert the dispersion model parameters
+ * into dictionary parameters for python.
+ */
+#include <vector>
+#include "dispersion_visitor.hh"
+
+using namespace std;
+
+/**
+ * Weight point class to hold averaging points
+ */
+class WeightPoint {
+public:
+	/// Value of the weight point
+	double value;
+	/// Weight of the weight point
+	double weight;
+
+	WeightPoint();
+	WeightPoint(double, double);
+};
+
+/**
+ * Basic averaging model. The class instance will
+ * generate a flat distribution of weight points
+ * according to the number of points specified
+ * and the width of the distribution. The center
+ * of the distribution is specified by the
+ * Parameter object taken in as a parameter.
+ */
+class DispersionModel {
+public:
+	/// Number of points to average over
+	int npts;
+	/// Width of the distribution (step function)
+	double width;
+
+	DispersionModel();
+	/// Method that generates the weight points
+	virtual void operator()(void *, vector<WeightPoint>&);
+	virtual void set_weights(int, double*, double*);
+	virtual void accept_as_source(DispersionVisitor*, void*, void*);
+	virtual void accept_as_destination(DispersionVisitor*, void*, void*);
+};
+
+
+/**
+ * Gaussian dispersion model
+ */
+class GaussianDispersion: public DispersionModel {
+public:
+	/// Number of sigmas on each side of the mean
+	int nsigmas;
+
+	GaussianDispersion();
+	void operator()(void *, vector<WeightPoint>&);
+	void accept_as_source(DispersionVisitor*, void*, void*);
+	void accept_as_destination(DispersionVisitor*, void*, void*);
+};
+
+/**
+ * Dispersion model based on arrays provided by the user
+ */
+class ArrayDispersion: public DispersionModel {
+private:
+	/// Array of values
+	double* _values;
+	/// Array of weights
+	double* _weights;
+
+	/// Method to set the weight points from arrays
+	void set_weights(int, double*, double*);
+	void operator()(void *, vector<WeightPoint>&);
+	void accept_as_source(DispersionVisitor*, void*, void*);
+	void accept_as_destination(DispersionVisitor*, void*, void*);
+public:
+
+};
+
+/**
+ * Parameter class to hold information about a
+ * parameter.
+ */
+class Parameter {
+public:
+	/// Current value of the parameter
+	double value;
+	/// True if the parameter has a minimum bound
+	bool has_min;
+	/// True if the parameter has a maximum bound
+	bool has_max;
+	/// Minimum bound
+	double min;
+	/// Maximum bound
+	double max;
+	/// True if the parameter can be dispersed or averaged
+	bool has_dispersion;
+	/// Pointer to the dispersion model object for this parameter
+	DispersionModel* dispersion;
+
+	Parameter();
+	Parameter(double);
+	Parameter(double, bool);
+
+	/// Method to set a minimum value
+	void set_min(double);
+	/// Method to set a maximum value
+	void set_max(double);
+	/// Method to get weight points for this parameter
+	void get_weights(vector<WeightPoint>&);
+	/// Returns the value of the parameter
+	double operator()();
+	/// Sets the value of the parameter
+	double operator=(double);
+};
+#endif
Index: sansmodels/src/sans/models/c_models/run.py
===================================================================
--- sansmodels/src/sans/models/c_models/run.py	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/run.py	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,20 @@
+import sys, os, math
+from sans.models.CylinderModel import CylinderModel
+from sans.models.DisperseModel import DisperseModel
+
+os.system('gcc -c ../../../libigor/libCylinder.c -o libigor.o -I../c_extensions -I../../../libigor')
+os.system('gcc -c ../c_extensions/cylinder.c -o libcyl.o -I../c_extensions -I../../../libigor')
+os.system('g++ -mno-cygwin cylinder.cpp parameters.cpp  dispersion_visitor.cpp libigor.o libcyl.o -D__MODELS_STANDALONE__ -I../c_extensions -I../../../libigor')
+os.system('a.exe')
+
+
+cyl = CylinderModel()
+cyl.setParam('cyl_phi', 0.0)
+cyl.setParam('cyl_theta', 0.0)
+
+disp = DisperseModel(cyl, ['radius'], [5])
+disp.setParam('n_pts', 100)
+print "" 
+print "Cyl (python): ", cyl.run(0.001)
+print "Disp(pyhton): ", disp.run(0.001)
+
Index: sansmodels/src/sans/models/c_models/test_cyl.py
===================================================================
--- sansmodels/src/sans/models/c_models/test_cyl.py	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
+++ sansmodels/src/sans/models/c_models/test_cyl.py	(revision fca6936ab3704433fa2327bcd55367e8ea28de72)
@@ -0,0 +1,50 @@
+import sys, os, math
+from sans.models.CylinderModel import CylinderModel
+from sans.models.DisperseModel import DisperseModel
+
+def old_cyl(q=0.001):
+
+    cyl = CylinderModel()
+    cyl.setParam('cyl_phi', 0.0)
+    cyl.setParam('cyl_theta', 0.0)
+    
+    disp = DisperseModel(cyl, ['radius'], [5])
+    disp.setParam('n_pts', 100)
+    print "" 
+    print "Old Cyl : ", cyl.run(q)
+    print "Old Disp: ", disp.run(q)
+    
+    
+    
+#----------------------------------------------------------------------
+    
+cyl = CylinderModel()
+print cyl.run(0.001)
+print cyl.run([0.001, .7854])
+print cyl.runXY([0.001,0.001])
+
+print ""
+print cyl.dispersion['disp_radius']
+
+cyl.dispersion['disp_radius']['width'] = 9.345
+
+q=0.001
+print "\n q=", q
+print "New Cyl : ", cyl.run(q)
+
+cyl.dispersion['disp_radius']['width']=5
+cyl.dispersion['disp_radius']['npts']=100
+print "New Disp: ", cyl.run(q)
+
+# Angular averaging
+values_ph = []
+values_th = []
+weights   = []
+for i in range(100):
+    values_ph.append(2.0*math.pi/99.0*i)
+    values_th.append(math.pi/99.0*i)
+    weights.append(1.0)
+
+
+
+old_cyl()
