Index: sansmodels/src/sans/models/c_models/CBinaryHSModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CBinaryHSModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CBinaryHSModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,12 +90,12 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"vol_frac_ls",Py_BuildValue("d",0.100000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.001000));
-        PyDict_SetItemString(self->params,"vol_frac_ss",Py_BuildValue("d",0.200000));
-        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"ls_sld",Py_BuildValue("d",0.000003));
-        PyDict_SetItemString(self->params,"ss_sld",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"s_radius",Py_BuildValue("d",25.000000));
-        PyDict_SetItemString(self->params,"l_radius",Py_BuildValue("d",100.000000));
+        PyDict_SetItemString(self->params,"vol_frac_ls",Py_BuildValue("d",0.10000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"vol_frac_ss",Py_BuildValue("d",0.20000000));
+        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.00000636));
+        PyDict_SetItemString(self->params,"ls_sld",Py_BuildValue("d",0.00000350));
+        PyDict_SetItemString(self->params,"ss_sld",Py_BuildValue("d",0.00000050));
+        PyDict_SetItemString(self->params,"s_radius",Py_BuildValue("d",25.00000000));
+        PyDict_SetItemString(self->params,"l_radius",Py_BuildValue("d",100.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CBinaryHSPSF11Model.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CBinaryHSPSF11Model.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CBinaryHSPSF11Model.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,12 +90,12 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"vol_frac_ls",Py_BuildValue("d",0.200000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.001000));
-        PyDict_SetItemString(self->params,"vol_frac_ss",Py_BuildValue("d",0.200000));
-        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"ls_sld",Py_BuildValue("d",0.000003));
-        PyDict_SetItemString(self->params,"ss_sld",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"s_radius",Py_BuildValue("d",25.000000));
-        PyDict_SetItemString(self->params,"l_radius",Py_BuildValue("d",160.000000));
+        PyDict_SetItemString(self->params,"vol_frac_ls",Py_BuildValue("d",0.20000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"vol_frac_ss",Py_BuildValue("d",0.20000000));
+        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.00000636));
+        PyDict_SetItemString(self->params,"ls_sld",Py_BuildValue("d",0.00000350));
+        PyDict_SetItemString(self->params,"ss_sld",Py_BuildValue("d",0.00000050));
+        PyDict_SetItemString(self->params,"s_radius",Py_BuildValue("d",25.00000000));
+        PyDict_SetItemString(self->params,"l_radius",Py_BuildValue("d",160.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CCoreFourShellModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CCoreFourShellModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/CCoreFourShellModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,645 @@
+/**
+	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
+ */
+
+/** CCoreFourShellModel
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY corefourshell.h
+ *          AND RE-RUN THE GENERATOR SCRIPT
+ *
+ */
+#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>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "corefourshell.h"
+}
+
+#include "models.hh"
+#include "dispersion_visitor.hh"
+
+/// Error object for raised exceptions
+static PyObject * CCoreFourShellModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Dispersion parameters
+    PyObject * dispersion;
+    /// Underlying model object
+    CoreFourShellModel * model;
+    /// Log for unit testing
+    PyObject * log;
+} CCoreFourShellModel;
+
+
+static void
+CCoreFourShellModel_dealloc(CCoreFourShellModel* self)
+{
+    Py_DECREF(self->params);
+    Py_DECREF(self->dispersion);
+    Py_DECREF(self->log);
+    delete self->model;
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CCoreFourShellModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CCoreFourShellModel *self;
+    
+    self = (CCoreFourShellModel *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CCoreFourShellModel_init(CCoreFourShellModel *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        self->dispersion = PyDict_New();
+        self->model = new CoreFourShellModel();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"thick_shell4",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"sld_shell4",Py_BuildValue("d",0.00000400));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sld_shell2",Py_BuildValue("d",0.00000200));
+        PyDict_SetItemString(self->params,"thick_shell1",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"thick_shell2",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"sld_shell1",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"rad_core",Py_BuildValue("d",60.00000000));
+        PyDict_SetItemString(self->params,"sld_core",Py_BuildValue("d",0.00000640));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"thick_shell3",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"sld_solv",Py_BuildValue("d",0.00000640));
+        PyDict_SetItemString(self->params,"sld_shell3",Py_BuildValue("d",0.00000300));
+        // Initialize dispersion / averaging parameter dict
+        DispersionVisitor* visitor = new DispersionVisitor();
+        PyObject * disp_dict;
+        disp_dict = PyDict_New();
+        self->model->rad_core.dispersion->accept_as_source(visitor, self->model->rad_core.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "rad_core", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->thick_shell1.dispersion->accept_as_source(visitor, self->model->thick_shell1.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "thick_shell1", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->thick_shell2.dispersion->accept_as_source(visitor, self->model->thick_shell2.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "thick_shell2", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->thick_shell3.dispersion->accept_as_source(visitor, self->model->thick_shell3.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "thick_shell3", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->thick_shell4.dispersion->accept_as_source(visitor, self->model->thick_shell4.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "thick_shell4", disp_dict);
+
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CCoreFourShellModel_members[] = {
+    {"params", T_OBJECT, offsetof(CCoreFourShellModel, params), 0,
+     "Parameters"},
+	{"dispersion", T_OBJECT, offsetof(CCoreFourShellModel, dispersion), 0,
+	  "Dispersion parameters"},     
+    {"log", T_OBJECT, offsetof(CCoreFourShellModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CCoreFourShellModel_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;
+    }
+}
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array q[] 
+ * @return: numpy array object 
+ */
+ 
+static PyObject *evaluateOneDim(CoreFourShellModel* 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  [x[],y[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDimXY( CoreFourShellModel* model, 
+                              PyArrayObject *x, PyArrayObject *y)
+ {
+    PyArrayObject *result;
+    int i,j, x_len, y_len, dims[1];
+    //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= "evaluateTwoDimXY  expect 2 numpy arrays";
+        PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+   
+	if (PyArray_Check(x) && PyArray_Check(y)) {
+		
+	    x_len = dims[0]= x->dimensions[0];
+        y_len = dims[0]= y->dimensions[0];
+	    
+	    // Make a new double matrix of same dims
+        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
+        if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	    }
+       
+        /* Do the calculation. */
+        for ( i=0; i< x_len; 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); 
+        
+        }else{
+		    PyErr_SetString(CCoreFourShellModelError, 
+                   "CCoreFourShellModel.evaluateTwoDimXY couldn't run.");
+	        return NULL;
+		}      	
+}
+/**
+ *  evalDistribution function evaluate a model function with input vector
+ *  @param args: input q as vector or [qx, qy] where qx, qy are vectors
+ *
+ */ 
+static PyObject * evalDistribution(CCoreFourShellModel *self, PyObject *args){
+	PyObject *qx, *qy;
+	PyArrayObject * pars;
+	int npars ,mpars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->thick_shell4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell4") );
+    self->model->sld_shell4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell4") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->sld_shell2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell2") );
+    self->model->thick_shell1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell1") );
+    self->model->thick_shell2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell2") );
+    self->model->sld_shell1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell1") );
+    self->model->rad_core = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rad_core") );
+    self->model->sld_core = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_core") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->thick_shell3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell3") );
+    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
+    self->model->sld_shell3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell3") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "rad_core");
+    self->model->rad_core.dispersion->accept_as_destination(visitor, self->model->rad_core.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell1");
+    self->model->thick_shell1.dispersion->accept_as_destination(visitor, self->model->thick_shell1.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell2");
+    self->model->thick_shell2.dispersion->accept_as_destination(visitor, self->model->thick_shell2.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell3");
+    self->model->thick_shell3.dispersion->accept_as_destination(visitor, self->model->thick_shell3.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell4");
+    self->model->thick_shell4.dispersion->accept_as_destination(visitor, self->model->thick_shell4.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(CCoreFourShellModelError, 
+	    	"CCoreFourShellModel.evalDistribution expects a q value.");
+		return NULL;
+	}
+    // Check params
+	
+    if(PyArray_Check(pars)==1) {
+		
+	    // Length of list should 1 or 2
+	    npars = pars->nd; 
+	    if(npars==1) {
+	        // input is a numpy array
+	        if (PyArray_Check(pars)) {
+		        return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+		    }
+		}else{
+		    PyErr_SetString(CCoreFourShellModelError, 
+                   "CCoreFourShellModel.evalDistribution expect numpy array of one dimension.");
+	        return NULL;
+		}
+    }else if( PyList_Check(pars)==1) {
+    	// Length of list should be 2 for I(qx,qy)
+	    mpars = PyList_GET_SIZE(pars); 
+	    if(mpars!=2) {
+	    	PyErr_SetString(CCoreFourShellModelError, 
+	    		"CCoreFourShellModel.evalDistribution expects a list of dimension 2.");
+	    	return NULL;
+	    }
+	     qx = PyList_GET_ITEM(pars,0);
+	     qy = PyList_GET_ITEM(pars,1);
+	     if (PyArray_Check(qx) && PyArray_Check(qy)) {
+	         return evaluateTwoDimXY(self->model, (PyArrayObject*)qx,
+		           (PyArrayObject*)qy);
+		 }else{
+		    PyErr_SetString(CCoreFourShellModelError, 
+                   "CCoreFourShellModel.evalDistribution expect 2 numpy arrays in list.");
+	        return NULL;
+	     }
+	}
+	PyErr_SetString(CCoreFourShellModelError, 
+                   "CCoreFourShellModel.evalDistribution couln't be run.");
+	return NULL;
+	
+}
+
+/**
+ * Function to call to evaluate model
+ * @param args: input q or [q,phi]
+ * @return: function value
+ */
+static PyObject * run(CCoreFourShellModel *self, PyObject *args) {
+	double q_value, phi_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->thick_shell4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell4") );
+    self->model->sld_shell4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell4") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->sld_shell2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell2") );
+    self->model->thick_shell1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell1") );
+    self->model->thick_shell2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell2") );
+    self->model->sld_shell1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell1") );
+    self->model->rad_core = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rad_core") );
+    self->model->sld_core = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_core") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->thick_shell3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell3") );
+    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
+    self->model->sld_shell3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell3") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "rad_core");
+    self->model->rad_core.dispersion->accept_as_destination(visitor, self->model->rad_core.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell1");
+    self->model->thick_shell1.dispersion->accept_as_destination(visitor, self->model->thick_shell1.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell2");
+    self->model->thick_shell2.dispersion->accept_as_destination(visitor, self->model->thick_shell2.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell3");
+    self->model->thick_shell3.dispersion->accept_as_destination(visitor, self->model->thick_shell3.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell4");
+    self->model->thick_shell4.dispersion->accept_as_destination(visitor, self->model->thick_shell4.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(CCoreFourShellModelError, 
+	    	"CCoreFourShellModel.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(CCoreFourShellModelError, 
+	    		"CCoreFourShellModel.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 = CCoreFourShellModel_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CCoreFourShellModel_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 = CCoreFourShellModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(q_value));
+	}	
+}
+/**
+ * Function to call to calculate_ER
+ * @return: effective radius value 
+ */
+static PyObject * calculate_ER(CCoreFourShellModel *self) {
+
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->thick_shell4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell4") );
+    self->model->sld_shell4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell4") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->sld_shell2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell2") );
+    self->model->thick_shell1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell1") );
+    self->model->thick_shell2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell2") );
+    self->model->sld_shell1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell1") );
+    self->model->rad_core = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rad_core") );
+    self->model->sld_core = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_core") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->thick_shell3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell3") );
+    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
+    self->model->sld_shell3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell3") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "rad_core");
+    self->model->rad_core.dispersion->accept_as_destination(visitor, self->model->rad_core.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell1");
+    self->model->thick_shell1.dispersion->accept_as_destination(visitor, self->model->thick_shell1.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell2");
+    self->model->thick_shell2.dispersion->accept_as_destination(visitor, self->model->thick_shell2.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell3");
+    self->model->thick_shell3.dispersion->accept_as_destination(visitor, self->model->thick_shell3.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell4");
+    self->model->thick_shell4.dispersion->accept_as_destination(visitor, self->model->thick_shell4.dispersion, disp_dict);
+
+		
+	return Py_BuildValue("d",(*(self->model)).calculate_ER());
+
+}
+/**
+ * Function to call to evaluate model in cartesian coordinates
+ * @param args: input q or [qx, qy]]
+ * @return: function value
+ */
+static PyObject * runXY(CCoreFourShellModel *self, PyObject *args) {
+	double qx_value, qy_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->thick_shell4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell4") );
+    self->model->sld_shell4 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell4") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->sld_shell2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell2") );
+    self->model->thick_shell1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell1") );
+    self->model->thick_shell2 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell2") );
+    self->model->sld_shell1 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell1") );
+    self->model->rad_core = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rad_core") );
+    self->model->sld_core = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_core") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->thick_shell3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thick_shell3") );
+    self->model->sld_solv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_solv") );
+    self->model->sld_shell3 = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_shell3") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "rad_core");
+    self->model->rad_core.dispersion->accept_as_destination(visitor, self->model->rad_core.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell1");
+    self->model->thick_shell1.dispersion->accept_as_destination(visitor, self->model->thick_shell1.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell2");
+    self->model->thick_shell2.dispersion->accept_as_destination(visitor, self->model->thick_shell2.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell3");
+    self->model->thick_shell3.dispersion->accept_as_destination(visitor, self->model->thick_shell3.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "thick_shell4");
+    self->model->thick_shell4.dispersion->accept_as_destination(visitor, self->model->thick_shell4.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(CCoreFourShellModelError, 
+	    	"CCoreFourShellModel.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(CCoreFourShellModelError, 
+	    		"CCoreFourShellModel.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 = CCoreFourShellModel_readDouble(PyList_GET_ITEM(pars,0));
+	    qy_value = CCoreFourShellModel_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 = CCoreFourShellModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(qx_value));
+	}	
+}
+
+static PyObject * reset(CCoreFourShellModel *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+static PyObject * set_dispersion(CCoreFourShellModel *self, PyObject *args) {
+	PyObject * disp;
+	const char * par_name;
+
+	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
+	    PyErr_SetString(CCoreFourShellModelError,
+	    	"CCoreFourShellModel.set_dispersion expects a DispersionModel object.");
+		return NULL;
+	}
+	void *temp = PyCObject_AsVoidPtr(disp);
+	DispersionModel * dispersion = static_cast<DispersionModel *>(temp);
+
+
+	// Ugliness necessary to go from python to C
+	    // TODO: refactor this
+    if (!strcmp(par_name, "rad_core")) {
+        self->model->rad_core.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "thick_shell1")) {
+        self->model->thick_shell1.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "thick_shell2")) {
+        self->model->thick_shell2.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "thick_shell3")) {
+        self->model->thick_shell3.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "thick_shell4")) {
+        self->model->thick_shell4.dispersion = dispersion;
+    } else {
+	    PyErr_SetString(CCoreFourShellModelError,
+	    	"CCoreFourShellModel.set_dispersion expects a valid parameter name.");
+		return NULL;
+	}
+
+	DispersionVisitor* visitor = new DispersionVisitor();
+	PyObject * disp_dict = PyDict_New();
+	dispersion->accept_as_source(visitor, dispersion, disp_dict);
+	PyDict_SetItemString(self->dispersion, par_name, disp_dict);
+    return Py_BuildValue("i",1);
+}
+
+
+static PyMethodDef CCoreFourShellModel_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"},
+    {"calculate_ER",      (PyCFunction)calculate_ER     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+      
+    {"evalDistribution",  (PyCFunction)evalDistribution , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy vector "},
+    {"reset",    (PyCFunction)reset   , METH_VARARGS,
+      "Reset pair correlation"},
+    {"set_dispersion",      (PyCFunction)set_dispersion     , METH_VARARGS,
+      "Set the dispersion model for a given parameter"},
+   {NULL}
+};
+
+static PyTypeObject CCoreFourShellModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CCoreFourShellModel",             /*tp_name*/
+    sizeof(CCoreFourShellModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CCoreFourShellModel_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*/
+    "CCoreFourShellModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CCoreFourShellModel_methods,             /* tp_methods */
+    CCoreFourShellModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CCoreFourShellModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CCoreFourShellModel_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 addCCoreFourShellModel(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CCoreFourShellModelType) < 0)
+        return;
+
+    Py_INCREF(&CCoreFourShellModelType);
+    PyModule_AddObject(module, "CCoreFourShellModel", (PyObject *)&CCoreFourShellModelType);
+    
+    d = PyModule_GetDict(module);
+    CCoreFourShellModelError = PyErr_NewException("CCoreFourShellModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CCoreFourShellModelError", CCoreFourShellModelError);
+}
+
Index: sansmodels/src/sans/models/c_models/CCoreShellCylinderModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CCoreShellCylinderModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CCoreShellCylinderModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,14 +90,14 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.570000));
-        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"thickness",Py_BuildValue("d",10.000000));
-        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"shell_sld",Py_BuildValue("d",0.000004));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.57000000));
+        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"thickness",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"shell_sld",Py_BuildValue("d",0.00000400));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CCoreShellEllipsoidModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CCoreShellEllipsoidModel.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/CCoreShellEllipsoidModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,15 +90,15 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sld_core",Py_BuildValue("d",0.000002));
-        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"polar_shell",Py_BuildValue("d",30.000000));
-        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"equat_shell",Py_BuildValue("d",250.000000));
-        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.001000));
-        PyDict_SetItemString(self->params,"equat_core",Py_BuildValue("d",200.000000));
-        PyDict_SetItemString(self->params,"polar_core",Py_BuildValue("d",20.000000));
-        PyDict_SetItemString(self->params,"sld_shell",Py_BuildValue("d",0.000001));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sld_core",Py_BuildValue("d",0.00000200));
+        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"polar_shell",Py_BuildValue("d",30.00000000));
+        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.00000630));
+        PyDict_SetItemString(self->params,"equat_shell",Py_BuildValue("d",250.00000000));
+        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"equat_core",Py_BuildValue("d",200.00000000));
+        PyDict_SetItemString(self->params,"polar_core",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"sld_shell",Py_BuildValue("d",0.00000100));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CCoreShellModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CCoreShellModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CCoreShellModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,11 +90,11 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"thickness",Py_BuildValue("d",10.000000));
-        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.000003));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",60.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"shell_sld",Py_BuildValue("d",0.000002));
+        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"thickness",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.00000300));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",60.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"shell_sld",Py_BuildValue("d",0.00000200));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CCylinderModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CCylinderModel.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/CCylinderModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,12 +90,12 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sldCyl",Py_BuildValue("d",0.000004));
-        PyDict_SetItemString(self->params,"cyl_theta",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.000000));
-        PyDict_SetItemString(self->params,"cyl_phi",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sldCyl",Py_BuildValue("d",0.00000400));
+        PyDict_SetItemString(self->params,"cyl_theta",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"cyl_phi",Py_BuildValue("d",1.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CDiamCylFunc.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CDiamCylFunc.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CDiamCylFunc.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,6 +90,6 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.000000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CDiamEllipFunc.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CDiamEllipFunc.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CDiamEllipFunc.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,6 +90,6 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"radius_b",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"radius_a",Py_BuildValue("d",20.000000));
+        PyDict_SetItemString(self->params,"radius_b",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"radius_a",Py_BuildValue("d",20.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CEllipsoidModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CEllipsoidModel.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/CEllipsoidModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,12 +90,12 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.570000));
-        PyDict_SetItemString(self->params,"radius_b",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"radius_a",Py_BuildValue("d",20.000000));
-        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"sldEll",Py_BuildValue("d",0.000004));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.57000000));
+        PyDict_SetItemString(self->params,"radius_b",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"radius_a",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"sldEll",Py_BuildValue("d",0.00000400));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CEllipticalCylinderModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CEllipticalCylinderModel.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/CEllipticalCylinderModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,14 +90,14 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sldCyl",Py_BuildValue("d",0.000004));
-        PyDict_SetItemString(self->params,"cyl_psi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"r_minor",Py_BuildValue("d",20.000000));
-        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"cyl_theta",Py_BuildValue("d",1.570000));
-        PyDict_SetItemString(self->params,"r_ratio",Py_BuildValue("d",1.500000));
-        PyDict_SetItemString(self->params,"cyl_phi",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sldCyl",Py_BuildValue("d",0.00000400));
+        PyDict_SetItemString(self->params,"cyl_psi",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"r_minor",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"cyl_theta",Py_BuildValue("d",1.57000000));
+        PyDict_SetItemString(self->params,"r_ratio",Py_BuildValue("d",1.50000000));
+        PyDict_SetItemString(self->params,"cyl_phi",Py_BuildValue("d",0.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CFlexCylEllipXModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CFlexCylEllipXModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/CFlexCylEllipXModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,607 @@
+/**
+	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
+ */
+
+/** CFlexCylEllipXModel
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY flexcyl_ellipX.h
+ *          AND RE-RUN THE GENERATOR SCRIPT
+ *
+ */
+#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>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "flexcyl_ellipX.h"
+}
+
+#include "models.hh"
+#include "dispersion_visitor.hh"
+
+/// Error object for raised exceptions
+static PyObject * CFlexCylEllipXModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Dispersion parameters
+    PyObject * dispersion;
+    /// Underlying model object
+    FlexCylEllipXModel * model;
+    /// Log for unit testing
+    PyObject * log;
+} CFlexCylEllipXModel;
+
+
+static void
+CFlexCylEllipXModel_dealloc(CFlexCylEllipXModel* self)
+{
+    Py_DECREF(self->params);
+    Py_DECREF(self->dispersion);
+    Py_DECREF(self->log);
+    delete self->model;
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CFlexCylEllipXModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CFlexCylEllipXModel *self;
+    
+    self = (CFlexCylEllipXModel *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CFlexCylEllipXModel_init(CFlexCylEllipXModel *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        self->dispersion = PyDict_New();
+        self->model = new FlexCylEllipXModel();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sldCyl",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"axis_ratio",Py_BuildValue("d",1.50000000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",1000.00000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00010000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000630));
+        PyDict_SetItemString(self->params,"kuhn_length",Py_BuildValue("d",100.00000000));
+        // Initialize dispersion / averaging parameter dict
+        DispersionVisitor* visitor = new DispersionVisitor();
+        PyObject * disp_dict;
+        disp_dict = PyDict_New();
+        self->model->length.dispersion->accept_as_source(visitor, self->model->length.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "length", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->kuhn_length.dispersion->accept_as_source(visitor, self->model->kuhn_length.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "kuhn_length", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->radius.dispersion->accept_as_source(visitor, self->model->radius.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "radius", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->axis_ratio.dispersion->accept_as_source(visitor, self->model->axis_ratio.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "axis_ratio", disp_dict);
+
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CFlexCylEllipXModel_members[] = {
+    {"params", T_OBJECT, offsetof(CFlexCylEllipXModel, params), 0,
+     "Parameters"},
+	{"dispersion", T_OBJECT, offsetof(CFlexCylEllipXModel, dispersion), 0,
+	  "Dispersion parameters"},     
+    {"log", T_OBJECT, offsetof(CFlexCylEllipXModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CFlexCylEllipXModel_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;
+    }
+}
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array q[] 
+ * @return: numpy array object 
+ */
+ 
+static PyObject *evaluateOneDim(FlexCylEllipXModel* 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  [x[],y[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDimXY( FlexCylEllipXModel* model, 
+                              PyArrayObject *x, PyArrayObject *y)
+ {
+    PyArrayObject *result;
+    int i,j, x_len, y_len, dims[1];
+    //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= "evaluateTwoDimXY  expect 2 numpy arrays";
+        PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+   
+	if (PyArray_Check(x) && PyArray_Check(y)) {
+		
+	    x_len = dims[0]= x->dimensions[0];
+        y_len = dims[0]= y->dimensions[0];
+	    
+	    // Make a new double matrix of same dims
+        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
+        if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	    }
+       
+        /* Do the calculation. */
+        for ( i=0; i< x_len; 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); 
+        
+        }else{
+		    PyErr_SetString(CFlexCylEllipXModelError, 
+                   "CFlexCylEllipXModel.evaluateTwoDimXY couldn't run.");
+	        return NULL;
+		}      	
+}
+/**
+ *  evalDistribution function evaluate a model function with input vector
+ *  @param args: input q as vector or [qx, qy] where qx, qy are vectors
+ *
+ */ 
+static PyObject * evalDistribution(CFlexCylEllipXModel *self, PyObject *args){
+	PyObject *qx, *qy;
+	PyArrayObject * pars;
+	int npars ,mpars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->sldCyl = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldCyl") );
+    self->model->axis_ratio = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_ratio") );
+    self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->kuhn_length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "kuhn_length") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "length");
+    self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "kuhn_length");
+    self->model->kuhn_length.dispersion->accept_as_destination(visitor, self->model->kuhn_length.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "axis_ratio");
+    self->model->axis_ratio.dispersion->accept_as_destination(visitor, self->model->axis_ratio.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(CFlexCylEllipXModelError, 
+	    	"CFlexCylEllipXModel.evalDistribution expects a q value.");
+		return NULL;
+	}
+    // Check params
+	
+    if(PyArray_Check(pars)==1) {
+		
+	    // Length of list should 1 or 2
+	    npars = pars->nd; 
+	    if(npars==1) {
+	        // input is a numpy array
+	        if (PyArray_Check(pars)) {
+		        return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+		    }
+		}else{
+		    PyErr_SetString(CFlexCylEllipXModelError, 
+                   "CFlexCylEllipXModel.evalDistribution expect numpy array of one dimension.");
+	        return NULL;
+		}
+    }else if( PyList_Check(pars)==1) {
+    	// Length of list should be 2 for I(qx,qy)
+	    mpars = PyList_GET_SIZE(pars); 
+	    if(mpars!=2) {
+	    	PyErr_SetString(CFlexCylEllipXModelError, 
+	    		"CFlexCylEllipXModel.evalDistribution expects a list of dimension 2.");
+	    	return NULL;
+	    }
+	     qx = PyList_GET_ITEM(pars,0);
+	     qy = PyList_GET_ITEM(pars,1);
+	     if (PyArray_Check(qx) && PyArray_Check(qy)) {
+	         return evaluateTwoDimXY(self->model, (PyArrayObject*)qx,
+		           (PyArrayObject*)qy);
+		 }else{
+		    PyErr_SetString(CFlexCylEllipXModelError, 
+                   "CFlexCylEllipXModel.evalDistribution expect 2 numpy arrays in list.");
+	        return NULL;
+	     }
+	}
+	PyErr_SetString(CFlexCylEllipXModelError, 
+                   "CFlexCylEllipXModel.evalDistribution couln't be run.");
+	return NULL;
+	
+}
+
+/**
+ * Function to call to evaluate model
+ * @param args: input q or [q,phi]
+ * @return: function value
+ */
+static PyObject * run(CFlexCylEllipXModel *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->sldCyl = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldCyl") );
+    self->model->axis_ratio = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_ratio") );
+    self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->kuhn_length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "kuhn_length") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "length");
+    self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "kuhn_length");
+    self->model->kuhn_length.dispersion->accept_as_destination(visitor, self->model->kuhn_length.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "axis_ratio");
+    self->model->axis_ratio.dispersion->accept_as_destination(visitor, self->model->axis_ratio.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(CFlexCylEllipXModelError, 
+	    	"CFlexCylEllipXModel.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(CFlexCylEllipXModelError, 
+	    		"CFlexCylEllipXModel.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 = CFlexCylEllipXModel_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CFlexCylEllipXModel_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 = CFlexCylEllipXModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(q_value));
+	}	
+}
+/**
+ * Function to call to calculate_ER
+ * @return: effective radius value 
+ */
+static PyObject * calculate_ER(CFlexCylEllipXModel *self) {
+
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->sldCyl = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldCyl") );
+    self->model->axis_ratio = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_ratio") );
+    self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->kuhn_length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "kuhn_length") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "length");
+    self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "kuhn_length");
+    self->model->kuhn_length.dispersion->accept_as_destination(visitor, self->model->kuhn_length.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "axis_ratio");
+    self->model->axis_ratio.dispersion->accept_as_destination(visitor, self->model->axis_ratio.dispersion, disp_dict);
+
+		
+	return Py_BuildValue("d",(*(self->model)).calculate_ER());
+
+}
+/**
+ * Function to call to evaluate model in cartesian coordinates
+ * @param args: input q or [qx, qy]]
+ * @return: function value
+ */
+static PyObject * runXY(CFlexCylEllipXModel *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->sldCyl = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldCyl") );
+    self->model->axis_ratio = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_ratio") );
+    self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->kuhn_length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "kuhn_length") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "length");
+    self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "kuhn_length");
+    self->model->kuhn_length.dispersion->accept_as_destination(visitor, self->model->kuhn_length.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "axis_ratio");
+    self->model->axis_ratio.dispersion->accept_as_destination(visitor, self->model->axis_ratio.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(CFlexCylEllipXModelError, 
+	    	"CFlexCylEllipXModel.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(CFlexCylEllipXModelError, 
+	    		"CFlexCylEllipXModel.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 = CFlexCylEllipXModel_readDouble(PyList_GET_ITEM(pars,0));
+	    qy_value = CFlexCylEllipXModel_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 = CFlexCylEllipXModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(qx_value));
+	}	
+}
+
+static PyObject * reset(CFlexCylEllipXModel *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+static PyObject * set_dispersion(CFlexCylEllipXModel *self, PyObject *args) {
+	PyObject * disp;
+	const char * par_name;
+
+	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
+	    PyErr_SetString(CFlexCylEllipXModelError,
+	    	"CFlexCylEllipXModel.set_dispersion expects a DispersionModel object.");
+		return NULL;
+	}
+	void *temp = PyCObject_AsVoidPtr(disp);
+	DispersionModel * dispersion = static_cast<DispersionModel *>(temp);
+
+
+	// Ugliness necessary to go from python to C
+	    // TODO: refactor this
+    if (!strcmp(par_name, "length")) {
+        self->model->length.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "kuhn_length")) {
+        self->model->kuhn_length.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "radius")) {
+        self->model->radius.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "axis_ratio")) {
+        self->model->axis_ratio.dispersion = dispersion;
+    } else {
+	    PyErr_SetString(CFlexCylEllipXModelError,
+	    	"CFlexCylEllipXModel.set_dispersion expects a valid parameter name.");
+		return NULL;
+	}
+
+	DispersionVisitor* visitor = new DispersionVisitor();
+	PyObject * disp_dict = PyDict_New();
+	dispersion->accept_as_source(visitor, dispersion, disp_dict);
+	PyDict_SetItemString(self->dispersion, par_name, disp_dict);
+    return Py_BuildValue("i",1);
+}
+
+
+static PyMethodDef CFlexCylEllipXModel_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"},
+    {"calculate_ER",      (PyCFunction)calculate_ER     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+      
+    {"evalDistribution",  (PyCFunction)evalDistribution , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy vector "},
+    {"reset",    (PyCFunction)reset   , METH_VARARGS,
+      "Reset pair correlation"},
+    {"set_dispersion",      (PyCFunction)set_dispersion     , METH_VARARGS,
+      "Set the dispersion model for a given parameter"},
+   {NULL}
+};
+
+static PyTypeObject CFlexCylEllipXModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CFlexCylEllipXModel",             /*tp_name*/
+    sizeof(CFlexCylEllipXModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CFlexCylEllipXModel_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*/
+    "CFlexCylEllipXModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CFlexCylEllipXModel_methods,             /* tp_methods */
+    CFlexCylEllipXModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CFlexCylEllipXModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CFlexCylEllipXModel_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 addCFlexCylEllipXModel(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CFlexCylEllipXModelType) < 0)
+        return;
+
+    Py_INCREF(&CFlexCylEllipXModelType);
+    PyModule_AddObject(module, "CFlexCylEllipXModel", (PyObject *)&CFlexCylEllipXModelType);
+    
+    d = PyModule_GetDict(module);
+    CFlexCylEllipXModelError = PyErr_NewException("CFlexCylEllipXModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CFlexCylEllipXModelError", CFlexCylEllipXModelError);
+}
+
Index: sansmodels/src/sans/models/c_models/CFlexibleCylinderModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CFlexibleCylinderModel.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/CFlexibleCylinderModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,11 +90,11 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sldCyl",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",1000.000000));
-        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000100));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.000000));
-        PyDict_SetItemString(self->params,"kuhn_length",Py_BuildValue("d",100.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sldCyl",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",1000.00000000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000630));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00010000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"kuhn_length",Py_BuildValue("d",100.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CFractalModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CFractalModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/CFractalModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,563 @@
+/**
+	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
+ */
+
+/** CFractalModel
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY fractal.h
+ *          AND RE-RUN THE GENERATOR SCRIPT
+ *
+ */
+#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>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "fractal.h"
+}
+
+#include "models.hh"
+#include "dispersion_visitor.hh"
+
+/// Error object for raised exceptions
+static PyObject * CFractalModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Dispersion parameters
+    PyObject * dispersion;
+    /// Underlying model object
+    FractalModel * model;
+    /// Log for unit testing
+    PyObject * log;
+} CFractalModel;
+
+
+static void
+CFractalModel_dealloc(CFractalModel* self)
+{
+    Py_DECREF(self->params);
+    Py_DECREF(self->dispersion);
+    Py_DECREF(self->log);
+    delete self->model;
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CFractalModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CFractalModel *self;
+    
+    self = (CFractalModel *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CFractalModel_init(CFractalModel *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        self->dispersion = PyDict_New();
+        self->model = new FractalModel();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",0.05000000));
+        PyDict_SetItemString(self->params,"fractal_dim",Py_BuildValue("d",2.00000000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000635));
+        PyDict_SetItemString(self->params,"cor_length",Py_BuildValue("d",100.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",5.00000000));
+        PyDict_SetItemString(self->params,"sldBlock",Py_BuildValue("d",0.00000200));
+        // Initialize dispersion / averaging parameter dict
+        DispersionVisitor* visitor = new DispersionVisitor();
+        PyObject * disp_dict;
+        disp_dict = PyDict_New();
+        self->model->radius.dispersion->accept_as_source(visitor, self->model->radius.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "radius", disp_dict);
+
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CFractalModel_members[] = {
+    {"params", T_OBJECT, offsetof(CFractalModel, params), 0,
+     "Parameters"},
+	{"dispersion", T_OBJECT, offsetof(CFractalModel, dispersion), 0,
+	  "Dispersion parameters"},     
+    {"log", T_OBJECT, offsetof(CFractalModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CFractalModel_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;
+    }
+}
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array q[] 
+ * @return: numpy array object 
+ */
+ 
+static PyObject *evaluateOneDim(FractalModel* 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  [x[],y[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDimXY( FractalModel* model, 
+                              PyArrayObject *x, PyArrayObject *y)
+ {
+    PyArrayObject *result;
+    int i,j, x_len, y_len, dims[1];
+    //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= "evaluateTwoDimXY  expect 2 numpy arrays";
+        PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+   
+	if (PyArray_Check(x) && PyArray_Check(y)) {
+		
+	    x_len = dims[0]= x->dimensions[0];
+        y_len = dims[0]= y->dimensions[0];
+	    
+	    // Make a new double matrix of same dims
+        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
+        if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	    }
+       
+        /* Do the calculation. */
+        for ( i=0; i< x_len; 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); 
+        
+        }else{
+		    PyErr_SetString(CFractalModelError, 
+                   "CFractalModel.evaluateTwoDimXY couldn't run.");
+	        return NULL;
+		}      	
+}
+/**
+ *  evalDistribution function evaluate a model function with input vector
+ *  @param args: input q as vector or [qx, qy] where qx, qy are vectors
+ *
+ */ 
+static PyObject * evalDistribution(CFractalModel *self, PyObject *args){
+	PyObject *qx, *qy;
+	PyArrayObject * pars;
+	int npars ,mpars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->fractal_dim = PyFloat_AsDouble( PyDict_GetItemString(self->params, "fractal_dim") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->cor_length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cor_length") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->sldBlock = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldBlock") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.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(CFractalModelError, 
+	    	"CFractalModel.evalDistribution expects a q value.");
+		return NULL;
+	}
+    // Check params
+	
+    if(PyArray_Check(pars)==1) {
+		
+	    // Length of list should 1 or 2
+	    npars = pars->nd; 
+	    if(npars==1) {
+	        // input is a numpy array
+	        if (PyArray_Check(pars)) {
+		        return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+		    }
+		}else{
+		    PyErr_SetString(CFractalModelError, 
+                   "CFractalModel.evalDistribution expect numpy array of one dimension.");
+	        return NULL;
+		}
+    }else if( PyList_Check(pars)==1) {
+    	// Length of list should be 2 for I(qx,qy)
+	    mpars = PyList_GET_SIZE(pars); 
+	    if(mpars!=2) {
+	    	PyErr_SetString(CFractalModelError, 
+	    		"CFractalModel.evalDistribution expects a list of dimension 2.");
+	    	return NULL;
+	    }
+	     qx = PyList_GET_ITEM(pars,0);
+	     qy = PyList_GET_ITEM(pars,1);
+	     if (PyArray_Check(qx) && PyArray_Check(qy)) {
+	         return evaluateTwoDimXY(self->model, (PyArrayObject*)qx,
+		           (PyArrayObject*)qy);
+		 }else{
+		    PyErr_SetString(CFractalModelError, 
+                   "CFractalModel.evalDistribution expect 2 numpy arrays in list.");
+	        return NULL;
+	     }
+	}
+	PyErr_SetString(CFractalModelError, 
+                   "CFractalModel.evalDistribution couln't be run.");
+	return NULL;
+	
+}
+
+/**
+ * Function to call to evaluate model
+ * @param args: input q or [q,phi]
+ * @return: function value
+ */
+static PyObject * run(CFractalModel *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->fractal_dim = PyFloat_AsDouble( PyDict_GetItemString(self->params, "fractal_dim") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->cor_length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cor_length") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->sldBlock = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldBlock") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.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(CFractalModelError, 
+	    	"CFractalModel.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(CFractalModelError, 
+	    		"CFractalModel.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 = CFractalModel_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CFractalModel_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 = CFractalModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(q_value));
+	}	
+}
+/**
+ * Function to call to calculate_ER
+ * @return: effective radius value 
+ */
+static PyObject * calculate_ER(CFractalModel *self) {
+
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->fractal_dim = PyFloat_AsDouble( PyDict_GetItemString(self->params, "fractal_dim") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->cor_length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cor_length") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->sldBlock = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldBlock") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+
+		
+	return Py_BuildValue("d",(*(self->model)).calculate_ER());
+
+}
+/**
+ * Function to call to evaluate model in cartesian coordinates
+ * @param args: input q or [qx, qy]]
+ * @return: function value
+ */
+static PyObject * runXY(CFractalModel *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->fractal_dim = PyFloat_AsDouble( PyDict_GetItemString(self->params, "fractal_dim") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->cor_length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cor_length") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->sldBlock = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldBlock") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.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(CFractalModelError, 
+	    	"CFractalModel.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(CFractalModelError, 
+	    		"CFractalModel.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 = CFractalModel_readDouble(PyList_GET_ITEM(pars,0));
+	    qy_value = CFractalModel_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 = CFractalModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(qx_value));
+	}	
+}
+
+static PyObject * reset(CFractalModel *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+static PyObject * set_dispersion(CFractalModel *self, PyObject *args) {
+	PyObject * disp;
+	const char * par_name;
+
+	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
+	    PyErr_SetString(CFractalModelError,
+	    	"CFractalModel.set_dispersion expects a DispersionModel object.");
+		return NULL;
+	}
+	void *temp = PyCObject_AsVoidPtr(disp);
+	DispersionModel * dispersion = static_cast<DispersionModel *>(temp);
+
+
+	// Ugliness necessary to go from python to C
+	    // TODO: refactor this
+    if (!strcmp(par_name, "radius")) {
+        self->model->radius.dispersion = dispersion;
+    } else {
+	    PyErr_SetString(CFractalModelError,
+	    	"CFractalModel.set_dispersion expects a valid parameter name.");
+		return NULL;
+	}
+
+	DispersionVisitor* visitor = new DispersionVisitor();
+	PyObject * disp_dict = PyDict_New();
+	dispersion->accept_as_source(visitor, dispersion, disp_dict);
+	PyDict_SetItemString(self->dispersion, par_name, disp_dict);
+    return Py_BuildValue("i",1);
+}
+
+
+static PyMethodDef CFractalModel_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"},
+    {"calculate_ER",      (PyCFunction)calculate_ER     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+      
+    {"evalDistribution",  (PyCFunction)evalDistribution , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy vector "},
+    {"reset",    (PyCFunction)reset   , METH_VARARGS,
+      "Reset pair correlation"},
+    {"set_dispersion",      (PyCFunction)set_dispersion     , METH_VARARGS,
+      "Set the dispersion model for a given parameter"},
+   {NULL}
+};
+
+static PyTypeObject CFractalModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CFractalModel",             /*tp_name*/
+    sizeof(CFractalModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CFractalModel_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*/
+    "CFractalModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CFractalModel_methods,             /* tp_methods */
+    CFractalModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CFractalModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CFractalModel_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 addCFractalModel(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CFractalModelType) < 0)
+        return;
+
+    Py_INCREF(&CFractalModelType);
+    PyModule_AddObject(module, "CFractalModel", (PyObject *)&CFractalModelType);
+    
+    d = PyModule_GetDict(module);
+    CFractalModelError = PyErr_NewException("CFractalModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CFractalModelError", CFractalModelError);
+}
+
Index: sansmodels/src/sans/models/c_models/CFuzzySphereModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CFuzzySphereModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/CFuzzySphereModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,571 @@
+/**
+	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
+ */
+
+/** CFuzzySphereModel
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY fuzzysphere.h
+ *          AND RE-RUN THE GENERATOR SCRIPT
+ *
+ */
+#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>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "fuzzysphere.h"
+}
+
+#include "models.hh"
+#include "dispersion_visitor.hh"
+
+/// Error object for raised exceptions
+static PyObject * CFuzzySphereModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Dispersion parameters
+    PyObject * dispersion;
+    /// Underlying model object
+    FuzzySphereModel * model;
+    /// Log for unit testing
+    PyObject * log;
+} CFuzzySphereModel;
+
+
+static void
+CFuzzySphereModel_dealloc(CFuzzySphereModel* self)
+{
+    Py_DECREF(self->params);
+    Py_DECREF(self->dispersion);
+    Py_DECREF(self->log);
+    delete self->model;
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CFuzzySphereModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CFuzzySphereModel *self;
+    
+    self = (CFuzzySphereModel *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CFuzzySphereModel_init(CFuzzySphereModel *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        self->dispersion = PyDict_New();
+        self->model = new FuzzySphereModel();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",0.01000000));
+        PyDict_SetItemString(self->params,"sldSph",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",60.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000300));
+        PyDict_SetItemString(self->params,"fuzziness",Py_BuildValue("d",10.00000000));
+        // Initialize dispersion / averaging parameter dict
+        DispersionVisitor* visitor = new DispersionVisitor();
+        PyObject * disp_dict;
+        disp_dict = PyDict_New();
+        self->model->radius.dispersion->accept_as_source(visitor, self->model->radius.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "radius", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->fuzziness.dispersion->accept_as_source(visitor, self->model->fuzziness.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "fuzziness", disp_dict);
+
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CFuzzySphereModel_members[] = {
+    {"params", T_OBJECT, offsetof(CFuzzySphereModel, params), 0,
+     "Parameters"},
+	{"dispersion", T_OBJECT, offsetof(CFuzzySphereModel, dispersion), 0,
+	  "Dispersion parameters"},     
+    {"log", T_OBJECT, offsetof(CFuzzySphereModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CFuzzySphereModel_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;
+    }
+}
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array q[] 
+ * @return: numpy array object 
+ */
+ 
+static PyObject *evaluateOneDim(FuzzySphereModel* 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  [x[],y[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDimXY( FuzzySphereModel* model, 
+                              PyArrayObject *x, PyArrayObject *y)
+ {
+    PyArrayObject *result;
+    int i,j, x_len, y_len, dims[1];
+    //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= "evaluateTwoDimXY  expect 2 numpy arrays";
+        PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+   
+	if (PyArray_Check(x) && PyArray_Check(y)) {
+		
+	    x_len = dims[0]= x->dimensions[0];
+        y_len = dims[0]= y->dimensions[0];
+	    
+	    // Make a new double matrix of same dims
+        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
+        if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	    }
+       
+        /* Do the calculation. */
+        for ( i=0; i< x_len; 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); 
+        
+        }else{
+		    PyErr_SetString(CFuzzySphereModelError, 
+                   "CFuzzySphereModel.evaluateTwoDimXY couldn't run.");
+	        return NULL;
+		}      	
+}
+/**
+ *  evalDistribution function evaluate a model function with input vector
+ *  @param args: input q as vector or [qx, qy] where qx, qy are vectors
+ *
+ */ 
+static PyObject * evalDistribution(CFuzzySphereModel *self, PyObject *args){
+	PyObject *qx, *qy;
+	PyArrayObject * pars;
+	int npars ,mpars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->sldSph = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSph") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->fuzziness = PyFloat_AsDouble( PyDict_GetItemString(self->params, "fuzziness") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "fuzziness");
+    self->model->fuzziness.dispersion->accept_as_destination(visitor, self->model->fuzziness.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(CFuzzySphereModelError, 
+	    	"CFuzzySphereModel.evalDistribution expects a q value.");
+		return NULL;
+	}
+    // Check params
+	
+    if(PyArray_Check(pars)==1) {
+		
+	    // Length of list should 1 or 2
+	    npars = pars->nd; 
+	    if(npars==1) {
+	        // input is a numpy array
+	        if (PyArray_Check(pars)) {
+		        return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+		    }
+		}else{
+		    PyErr_SetString(CFuzzySphereModelError, 
+                   "CFuzzySphereModel.evalDistribution expect numpy array of one dimension.");
+	        return NULL;
+		}
+    }else if( PyList_Check(pars)==1) {
+    	// Length of list should be 2 for I(qx,qy)
+	    mpars = PyList_GET_SIZE(pars); 
+	    if(mpars!=2) {
+	    	PyErr_SetString(CFuzzySphereModelError, 
+	    		"CFuzzySphereModel.evalDistribution expects a list of dimension 2.");
+	    	return NULL;
+	    }
+	     qx = PyList_GET_ITEM(pars,0);
+	     qy = PyList_GET_ITEM(pars,1);
+	     if (PyArray_Check(qx) && PyArray_Check(qy)) {
+	         return evaluateTwoDimXY(self->model, (PyArrayObject*)qx,
+		           (PyArrayObject*)qy);
+		 }else{
+		    PyErr_SetString(CFuzzySphereModelError, 
+                   "CFuzzySphereModel.evalDistribution expect 2 numpy arrays in list.");
+	        return NULL;
+	     }
+	}
+	PyErr_SetString(CFuzzySphereModelError, 
+                   "CFuzzySphereModel.evalDistribution couln't be run.");
+	return NULL;
+	
+}
+
+/**
+ * Function to call to evaluate model
+ * @param args: input q or [q,phi]
+ * @return: function value
+ */
+static PyObject * run(CFuzzySphereModel *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->sldSph = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSph") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->fuzziness = PyFloat_AsDouble( PyDict_GetItemString(self->params, "fuzziness") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "fuzziness");
+    self->model->fuzziness.dispersion->accept_as_destination(visitor, self->model->fuzziness.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(CFuzzySphereModelError, 
+	    	"CFuzzySphereModel.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(CFuzzySphereModelError, 
+	    		"CFuzzySphereModel.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 = CFuzzySphereModel_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CFuzzySphereModel_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 = CFuzzySphereModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(q_value));
+	}	
+}
+/**
+ * Function to call to calculate_ER
+ * @return: effective radius value 
+ */
+static PyObject * calculate_ER(CFuzzySphereModel *self) {
+
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->sldSph = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSph") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->fuzziness = PyFloat_AsDouble( PyDict_GetItemString(self->params, "fuzziness") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "fuzziness");
+    self->model->fuzziness.dispersion->accept_as_destination(visitor, self->model->fuzziness.dispersion, disp_dict);
+
+		
+	return Py_BuildValue("d",(*(self->model)).calculate_ER());
+
+}
+/**
+ * Function to call to evaluate model in cartesian coordinates
+ * @param args: input q or [qx, qy]]
+ * @return: function value
+ */
+static PyObject * runXY(CFuzzySphereModel *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->sldSph = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSph") );
+    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
+    self->model->fuzziness = PyFloat_AsDouble( PyDict_GetItemString(self->params, "fuzziness") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
+    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "fuzziness");
+    self->model->fuzziness.dispersion->accept_as_destination(visitor, self->model->fuzziness.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(CFuzzySphereModelError, 
+	    	"CFuzzySphereModel.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(CFuzzySphereModelError, 
+	    		"CFuzzySphereModel.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 = CFuzzySphereModel_readDouble(PyList_GET_ITEM(pars,0));
+	    qy_value = CFuzzySphereModel_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 = CFuzzySphereModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(qx_value));
+	}	
+}
+
+static PyObject * reset(CFuzzySphereModel *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+static PyObject * set_dispersion(CFuzzySphereModel *self, PyObject *args) {
+	PyObject * disp;
+	const char * par_name;
+
+	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
+	    PyErr_SetString(CFuzzySphereModelError,
+	    	"CFuzzySphereModel.set_dispersion expects a DispersionModel object.");
+		return NULL;
+	}
+	void *temp = PyCObject_AsVoidPtr(disp);
+	DispersionModel * dispersion = static_cast<DispersionModel *>(temp);
+
+
+	// Ugliness necessary to go from python to C
+	    // TODO: refactor this
+    if (!strcmp(par_name, "radius")) {
+        self->model->radius.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "fuzziness")) {
+        self->model->fuzziness.dispersion = dispersion;
+    } else {
+	    PyErr_SetString(CFuzzySphereModelError,
+	    	"CFuzzySphereModel.set_dispersion expects a valid parameter name.");
+		return NULL;
+	}
+
+	DispersionVisitor* visitor = new DispersionVisitor();
+	PyObject * disp_dict = PyDict_New();
+	dispersion->accept_as_source(visitor, dispersion, disp_dict);
+	PyDict_SetItemString(self->dispersion, par_name, disp_dict);
+    return Py_BuildValue("i",1);
+}
+
+
+static PyMethodDef CFuzzySphereModel_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"},
+    {"calculate_ER",      (PyCFunction)calculate_ER     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+      
+    {"evalDistribution",  (PyCFunction)evalDistribution , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy vector "},
+    {"reset",    (PyCFunction)reset   , METH_VARARGS,
+      "Reset pair correlation"},
+    {"set_dispersion",      (PyCFunction)set_dispersion     , METH_VARARGS,
+      "Set the dispersion model for a given parameter"},
+   {NULL}
+};
+
+static PyTypeObject CFuzzySphereModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CFuzzySphereModel",             /*tp_name*/
+    sizeof(CFuzzySphereModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CFuzzySphereModel_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*/
+    "CFuzzySphereModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CFuzzySphereModel_methods,             /* tp_methods */
+    CFuzzySphereModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CFuzzySphereModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CFuzzySphereModel_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 addCFuzzySphereModel(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CFuzzySphereModelType) < 0)
+        return;
+
+    Py_INCREF(&CFuzzySphereModelType);
+    PyModule_AddObject(module, "CFuzzySphereModel", (PyObject *)&CFuzzySphereModelType);
+    
+    d = PyModule_GetDict(module);
+    CFuzzySphereModelError = PyErr_NewException("CFuzzySphereModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CFuzzySphereModelError", CFuzzySphereModelError);
+}
+
Index: sansmodels/src/sans/models/c_models/CGaussian.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CGaussian.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CGaussian.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,7 +90,7 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sigma",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"center",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sigma",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"center",Py_BuildValue("d",0.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CHardsphereStructure.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CHardsphereStructure.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CHardsphereStructure.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,6 +90,6 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"effect_radius",Py_BuildValue("d",50.000000));
-        PyDict_SetItemString(self->params,"volfraction",Py_BuildValue("d",0.200000));
+        PyDict_SetItemString(self->params,"effect_radius",Py_BuildValue("d",50.00000000));
+        PyDict_SetItemString(self->params,"volfraction",Py_BuildValue("d",0.20000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CHayterMSAStructure.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CHayterMSAStructure.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CHayterMSAStructure.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,10 +90,10 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"temperature",Py_BuildValue("d",318.160000));
-        PyDict_SetItemString(self->params,"saltconc",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"dielectconst",Py_BuildValue("d",71.080000));
-        PyDict_SetItemString(self->params,"volfraction",Py_BuildValue("d",0.019200));
-        PyDict_SetItemString(self->params,"charge",Py_BuildValue("d",19.000000));
-        PyDict_SetItemString(self->params,"effect_radius",Py_BuildValue("d",20.750000));
+        PyDict_SetItemString(self->params,"temperature",Py_BuildValue("d",318.16000000));
+        PyDict_SetItemString(self->params,"saltconc",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"dielectconst",Py_BuildValue("d",71.08000000));
+        PyDict_SetItemString(self->params,"volfraction",Py_BuildValue("d",0.01920000));
+        PyDict_SetItemString(self->params,"charge",Py_BuildValue("d",19.00000000));
+        PyDict_SetItemString(self->params,"effect_radius",Py_BuildValue("d",20.75000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CHollowCylinderModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CHollowCylinderModel.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/CHollowCylinderModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,13 +90,13 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sldCyl",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"core_radius",Py_BuildValue("d",20.000000));
-        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.570000));
-        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.010000));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",30.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sldCyl",Py_BuildValue("d",0.00000630));
+        PyDict_SetItemString(self->params,"core_radius",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.57000000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.01000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",30.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CLamellarFFHGModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CLamellarFFHGModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CLamellarFFHGModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,11 +90,11 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"h_thickness",Py_BuildValue("d",10.000000));
-        PyDict_SetItemString(self->params,"sld_tail",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"t_length",Py_BuildValue("d",15.000000));
-        PyDict_SetItemString(self->params,"sld_head",Py_BuildValue("d",0.000003));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"h_thickness",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"sld_tail",Py_BuildValue("d",0.00000040));
+        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.00000600));
+        PyDict_SetItemString(self->params,"t_length",Py_BuildValue("d",15.00000000));
+        PyDict_SetItemString(self->params,"sld_head",Py_BuildValue("d",0.00000300));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CLamellarModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CLamellarModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CLamellarModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,9 +90,9 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"sld_sol",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"bi_thick",Py_BuildValue("d",50.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"sld_bi",Py_BuildValue("d",0.000001));
+        PyDict_SetItemString(self->params,"sld_sol",Py_BuildValue("d",0.00000630));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"bi_thick",Py_BuildValue("d",50.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"sld_bi",Py_BuildValue("d",0.00000100));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CLamellarPSHGModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CLamellarPSHGModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CLamellarPSHGModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,14 +90,14 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"n_plates",Py_BuildValue("d",30.000000));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"deltaT",Py_BuildValue("d",10.000000));
-        PyDict_SetItemString(self->params,"spacing",Py_BuildValue("d",40.000000));
-        PyDict_SetItemString(self->params,"sld_tail",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"caille",Py_BuildValue("d",0.001000));
-        PyDict_SetItemString(self->params,"sld_head",Py_BuildValue("d",0.000002));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.001000));
-        PyDict_SetItemString(self->params,"deltaH",Py_BuildValue("d",2.000000));
+        PyDict_SetItemString(self->params,"n_plates",Py_BuildValue("d",30.00000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"deltaT",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"spacing",Py_BuildValue("d",40.00000000));
+        PyDict_SetItemString(self->params,"sld_tail",Py_BuildValue("d",0.00000040));
+        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.00000600));
+        PyDict_SetItemString(self->params,"caille",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"sld_head",Py_BuildValue("d",0.00000200));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"deltaH",Py_BuildValue("d",2.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CLamellarPSModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CLamellarPSModel.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/CLamellarPSModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,12 +90,12 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"n_plates",Py_BuildValue("d",20.000000));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sld_sol",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"delta",Py_BuildValue("d",30.000000));
-        PyDict_SetItemString(self->params,"spacing",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"caille",Py_BuildValue("d",0.100000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"sld_bi",Py_BuildValue("d",0.000006));
+        PyDict_SetItemString(self->params,"n_plates",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sld_sol",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"delta",Py_BuildValue("d",30.00000000));
+        PyDict_SetItemString(self->params,"spacing",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"caille",Py_BuildValue("d",0.10000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"sld_bi",Py_BuildValue("d",0.00000630));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CLogNormal.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CLogNormal.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CLogNormal.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,7 +90,7 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sigma",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"center",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sigma",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"center",Py_BuildValue("d",0.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CLorentzian.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CLorentzian.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CLorentzian.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,7 +90,7 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"center",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"gamma",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"center",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"gamma",Py_BuildValue("d",1.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CMultiShellModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CMultiShellModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CMultiShellModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,12 +90,12 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"core_radius",Py_BuildValue("d",60.000000));
-        PyDict_SetItemString(self->params,"n_pairs",Py_BuildValue("d",2.000000));
-        PyDict_SetItemString(self->params,"w_thickness",Py_BuildValue("d",10.000000));
-        PyDict_SetItemString(self->params,"s_thickness",Py_BuildValue("d",10.000000));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"shell_sld",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.00000640));
+        PyDict_SetItemString(self->params,"core_radius",Py_BuildValue("d",60.00000000));
+        PyDict_SetItemString(self->params,"n_pairs",Py_BuildValue("d",2.00000000));
+        PyDict_SetItemString(self->params,"w_thickness",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"s_thickness",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"shell_sld",Py_BuildValue("d",0.00000040));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/COblateModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/COblateModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/COblateModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,14 +90,14 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"major_core",Py_BuildValue("d",200.000000));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"minor_core",Py_BuildValue("d",20.000000));
-        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.001000));
-        PyDict_SetItemString(self->params,"major_shell",Py_BuildValue("d",250.000000));
-        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"minor_shell",Py_BuildValue("d",30.000000));
+        PyDict_SetItemString(self->params,"major_core",Py_BuildValue("d",200.00000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"minor_core",Py_BuildValue("d",20.00000000));
+        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.00000630));
+        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"major_shell",Py_BuildValue("d",250.00000000));
+        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"minor_shell",Py_BuildValue("d",30.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CParallelepipedModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CParallelepipedModel.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/CParallelepipedModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,14 +90,14 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"short_a",Py_BuildValue("d",35.000000));
-        PyDict_SetItemString(self->params,"short_b",Py_BuildValue("d",75.000000));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"long_c",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"sldPipe",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"parallel_psi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"parallel_phi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"parallel_theta",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"short_a",Py_BuildValue("d",35.00000000));
+        PyDict_SetItemString(self->params,"short_b",Py_BuildValue("d",75.00000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"long_c",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"sldPipe",Py_BuildValue("d",0.00000630));
+        PyDict_SetItemString(self->params,"parallel_psi",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"parallel_phi",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"parallel_theta",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CPoly_GaussCoil.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CPoly_GaussCoil.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/CPoly_GaussCoil.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,548 @@
+/**
+	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
+ */
+
+/** CPoly_GaussCoil
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY polygausscoil.h
+ *          AND RE-RUN THE GENERATOR SCRIPT
+ *
+ */
+#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>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "polygausscoil.h"
+}
+
+#include "models.hh"
+#include "dispersion_visitor.hh"
+
+/// Error object for raised exceptions
+static PyObject * CPoly_GaussCoilError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Dispersion parameters
+    PyObject * dispersion;
+    /// Underlying model object
+    Poly_GaussCoil * model;
+    /// Log for unit testing
+    PyObject * log;
+} CPoly_GaussCoil;
+
+
+static void
+CPoly_GaussCoil_dealloc(CPoly_GaussCoil* self)
+{
+    Py_DECREF(self->params);
+    Py_DECREF(self->dispersion);
+    Py_DECREF(self->log);
+    delete self->model;
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CPoly_GaussCoil_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CPoly_GaussCoil *self;
+    
+    self = (CPoly_GaussCoil *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CPoly_GaussCoil_init(CPoly_GaussCoil *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        self->dispersion = PyDict_New();
+        self->model = new Poly_GaussCoil();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"poly_m",Py_BuildValue("d",2.00000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"rg",Py_BuildValue("d",60.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        // Initialize dispersion / averaging parameter dict
+        DispersionVisitor* visitor = new DispersionVisitor();
+        PyObject * disp_dict;
+        disp_dict = PyDict_New();
+        self->model->rg.dispersion->accept_as_source(visitor, self->model->rg.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "rg", disp_dict);
+
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CPoly_GaussCoil_members[] = {
+    {"params", T_OBJECT, offsetof(CPoly_GaussCoil, params), 0,
+     "Parameters"},
+	{"dispersion", T_OBJECT, offsetof(CPoly_GaussCoil, dispersion), 0,
+	  "Dispersion parameters"},     
+    {"log", T_OBJECT, offsetof(CPoly_GaussCoil, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CPoly_GaussCoil_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;
+    }
+}
+/**
+ * Function to call to evaluate model
+ * @param args: input numpy array q[] 
+ * @return: numpy array object 
+ */
+ 
+static PyObject *evaluateOneDim(Poly_GaussCoil* 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  [x[],y[]]
+ * @return: numpy array object 
+ */
+ static PyObject * evaluateTwoDimXY( Poly_GaussCoil* model, 
+                              PyArrayObject *x, PyArrayObject *y)
+ {
+    PyArrayObject *result;
+    int i,j, x_len, y_len, dims[1];
+    //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= "evaluateTwoDimXY  expect 2 numpy arrays";
+        PyErr_SetString(PyExc_ValueError , message); 
+        return NULL;
+    }
+   
+	if (PyArray_Check(x) && PyArray_Check(y)) {
+		
+	    x_len = dims[0]= x->dimensions[0];
+        y_len = dims[0]= y->dimensions[0];
+	    
+	    // Make a new double matrix of same dims
+        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
+        if (result == NULL){
+	    const char * message= "Could not create result ";
+        PyErr_SetString(PyExc_RuntimeError , message);
+	    return NULL;
+	    }
+       
+        /* Do the calculation. */
+        for ( i=0; i< x_len; 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); 
+        
+        }else{
+		    PyErr_SetString(CPoly_GaussCoilError, 
+                   "CPoly_GaussCoil.evaluateTwoDimXY couldn't run.");
+	        return NULL;
+		}      	
+}
+/**
+ *  evalDistribution function evaluate a model function with input vector
+ *  @param args: input q as vector or [qx, qy] where qx, qy are vectors
+ *
+ */ 
+static PyObject * evalDistribution(CPoly_GaussCoil *self, PyObject *args){
+	PyObject *qx, *qy;
+	PyArrayObject * pars;
+	int npars ,mpars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->poly_m = PyFloat_AsDouble( PyDict_GetItemString(self->params, "poly_m") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->rg = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rg") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "rg");
+    self->model->rg.dispersion->accept_as_destination(visitor, self->model->rg.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(CPoly_GaussCoilError, 
+	    	"CPoly_GaussCoil.evalDistribution expects a q value.");
+		return NULL;
+	}
+    // Check params
+	
+    if(PyArray_Check(pars)==1) {
+		
+	    // Length of list should 1 or 2
+	    npars = pars->nd; 
+	    if(npars==1) {
+	        // input is a numpy array
+	        if (PyArray_Check(pars)) {
+		        return evaluateOneDim(self->model, (PyArrayObject*)pars); 
+		    }
+		}else{
+		    PyErr_SetString(CPoly_GaussCoilError, 
+                   "CPoly_GaussCoil.evalDistribution expect numpy array of one dimension.");
+	        return NULL;
+		}
+    }else if( PyList_Check(pars)==1) {
+    	// Length of list should be 2 for I(qx,qy)
+	    mpars = PyList_GET_SIZE(pars); 
+	    if(mpars!=2) {
+	    	PyErr_SetString(CPoly_GaussCoilError, 
+	    		"CPoly_GaussCoil.evalDistribution expects a list of dimension 2.");
+	    	return NULL;
+	    }
+	     qx = PyList_GET_ITEM(pars,0);
+	     qy = PyList_GET_ITEM(pars,1);
+	     if (PyArray_Check(qx) && PyArray_Check(qy)) {
+	         return evaluateTwoDimXY(self->model, (PyArrayObject*)qx,
+		           (PyArrayObject*)qy);
+		 }else{
+		    PyErr_SetString(CPoly_GaussCoilError, 
+                   "CPoly_GaussCoil.evalDistribution expect 2 numpy arrays in list.");
+	        return NULL;
+	     }
+	}
+	PyErr_SetString(CPoly_GaussCoilError, 
+                   "CPoly_GaussCoil.evalDistribution couln't be run.");
+	return NULL;
+	
+}
+
+/**
+ * Function to call to evaluate model
+ * @param args: input q or [q,phi]
+ * @return: function value
+ */
+static PyObject * run(CPoly_GaussCoil *self, PyObject *args) {
+	double q_value, phi_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->poly_m = PyFloat_AsDouble( PyDict_GetItemString(self->params, "poly_m") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->rg = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rg") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "rg");
+    self->model->rg.dispersion->accept_as_destination(visitor, self->model->rg.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(CPoly_GaussCoilError, 
+	    	"CPoly_GaussCoil.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(CPoly_GaussCoilError, 
+	    		"CPoly_GaussCoil.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 = CPoly_GaussCoil_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CPoly_GaussCoil_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 = CPoly_GaussCoil_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(q_value));
+	}	
+}
+/**
+ * Function to call to calculate_ER
+ * @return: effective radius value 
+ */
+static PyObject * calculate_ER(CPoly_GaussCoil *self) {
+
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->poly_m = PyFloat_AsDouble( PyDict_GetItemString(self->params, "poly_m") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->rg = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rg") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "rg");
+    self->model->rg.dispersion->accept_as_destination(visitor, self->model->rg.dispersion, disp_dict);
+
+		
+	return Py_BuildValue("d",(*(self->model)).calculate_ER());
+
+}
+/**
+ * Function to call to evaluate model in cartesian coordinates
+ * @param args: input q or [qx, qy]]
+ * @return: function value
+ */
+static PyObject * runXY(CPoly_GaussCoil *self, PyObject *args) {
+	double qx_value, qy_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	    // Reader parameter dictionary
+    self->model->poly_m = PyFloat_AsDouble( PyDict_GetItemString(self->params, "poly_m") );
+    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model->rg = PyFloat_AsDouble( PyDict_GetItemString(self->params, "rg") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "rg");
+    self->model->rg.dispersion->accept_as_destination(visitor, self->model->rg.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(CPoly_GaussCoilError, 
+	    	"CPoly_GaussCoil.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(CPoly_GaussCoilError, 
+	    		"CPoly_GaussCoil.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 = CPoly_GaussCoil_readDouble(PyList_GET_ITEM(pars,0));
+	    qy_value = CPoly_GaussCoil_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 = CPoly_GaussCoil_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(qx_value));
+	}	
+}
+
+static PyObject * reset(CPoly_GaussCoil *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+static PyObject * set_dispersion(CPoly_GaussCoil *self, PyObject *args) {
+	PyObject * disp;
+	const char * par_name;
+
+	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
+	    PyErr_SetString(CPoly_GaussCoilError,
+	    	"CPoly_GaussCoil.set_dispersion expects a DispersionModel object.");
+		return NULL;
+	}
+	void *temp = PyCObject_AsVoidPtr(disp);
+	DispersionModel * dispersion = static_cast<DispersionModel *>(temp);
+
+
+	// Ugliness necessary to go from python to C
+	    // TODO: refactor this
+    if (!strcmp(par_name, "rg")) {
+        self->model->rg.dispersion = dispersion;
+    } else {
+	    PyErr_SetString(CPoly_GaussCoilError,
+	    	"CPoly_GaussCoil.set_dispersion expects a valid parameter name.");
+		return NULL;
+	}
+
+	DispersionVisitor* visitor = new DispersionVisitor();
+	PyObject * disp_dict = PyDict_New();
+	dispersion->accept_as_source(visitor, dispersion, disp_dict);
+	PyDict_SetItemString(self->dispersion, par_name, disp_dict);
+    return Py_BuildValue("i",1);
+}
+
+
+static PyMethodDef CPoly_GaussCoil_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"},
+    {"calculate_ER",      (PyCFunction)calculate_ER     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+      
+    {"evalDistribution",  (PyCFunction)evalDistribution , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy vector "},
+    {"reset",    (PyCFunction)reset   , METH_VARARGS,
+      "Reset pair correlation"},
+    {"set_dispersion",      (PyCFunction)set_dispersion     , METH_VARARGS,
+      "Set the dispersion model for a given parameter"},
+   {NULL}
+};
+
+static PyTypeObject CPoly_GaussCoilType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CPoly_GaussCoil",             /*tp_name*/
+    sizeof(CPoly_GaussCoil),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CPoly_GaussCoil_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*/
+    "CPoly_GaussCoil objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CPoly_GaussCoil_methods,             /* tp_methods */
+    CPoly_GaussCoil_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CPoly_GaussCoil_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CPoly_GaussCoil_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 addCPoly_GaussCoil(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CPoly_GaussCoilType) < 0)
+        return;
+
+    Py_INCREF(&CPoly_GaussCoilType);
+    PyModule_AddObject(module, "CPoly_GaussCoil", (PyObject *)&CPoly_GaussCoilType);
+    
+    d = PyModule_GetDict(module);
+    CPoly_GaussCoilError = PyErr_NewException("CPoly_GaussCoil.error", NULL, NULL);
+    PyDict_SetItemString(d, "CPoly_GaussCoilError", CPoly_GaussCoilError);
+}
+
Index: sansmodels/src/sans/models/c_models/CProlateModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CProlateModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CProlateModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,12 +90,12 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"major_core",Py_BuildValue("d",100.000000));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"minor_core",Py_BuildValue("d",50.000000));
-        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.001000));
-        PyDict_SetItemString(self->params,"major_shell",Py_BuildValue("d",110.000000));
-        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"minor_shell",Py_BuildValue("d",60.000000));
+        PyDict_SetItemString(self->params,"major_core",Py_BuildValue("d",100.00000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"minor_core",Py_BuildValue("d",50.00000000));
+        PyDict_SetItemString(self->params,"sld_solvent",Py_BuildValue("d",0.00000630));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"major_shell",Py_BuildValue("d",110.00000000));
+        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"minor_shell",Py_BuildValue("d",60.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CSchulz.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CSchulz.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CSchulz.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,7 +90,7 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"sigma",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"center",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"sigma",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"center",Py_BuildValue("d",0.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CSphereModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CSphereModel.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/CSphereModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,9 +90,9 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.000001));
-        PyDict_SetItemString(self->params,"sldSph",Py_BuildValue("d",0.000002));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",60.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.00000100));
+        PyDict_SetItemString(self->params,"sldSph",Py_BuildValue("d",0.00000200));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",60.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CSquareWellStructure.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CSquareWellStructure.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CSquareWellStructure.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,8 +90,8 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"welldepth",Py_BuildValue("d",1.500000));
-        PyDict_SetItemString(self->params,"effect_radius",Py_BuildValue("d",50.000000));
-        PyDict_SetItemString(self->params,"volfraction",Py_BuildValue("d",0.040000));
-        PyDict_SetItemString(self->params,"wellwidth",Py_BuildValue("d",1.200000));
+        PyDict_SetItemString(self->params,"welldepth",Py_BuildValue("d",1.50000000));
+        PyDict_SetItemString(self->params,"effect_radius",Py_BuildValue("d",50.00000000));
+        PyDict_SetItemString(self->params,"volfraction",Py_BuildValue("d",0.04000000));
+        PyDict_SetItemString(self->params,"wellwidth",Py_BuildValue("d",1.20000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CStackedDisksModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CStackedDisksModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CStackedDisksModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,16 +90,16 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.000004));
-        PyDict_SetItemString(self->params,"core_thick",Py_BuildValue("d",10.000000));
-        PyDict_SetItemString(self->params,"layer_thick",Py_BuildValue("d",15.000000));
-        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"layer_sld",Py_BuildValue("d",-0.000000));
-        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.000005));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",0.010000));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",3000.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.001000));
-        PyDict_SetItemString(self->params,"sigma_d",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"n_stacking",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.00000400));
+        PyDict_SetItemString(self->params,"core_thick",Py_BuildValue("d",10.00000000));
+        PyDict_SetItemString(self->params,"layer_thick",Py_BuildValue("d",15.00000000));
+        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"layer_sld",Py_BuildValue("d",-0.00000040));
+        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.00000500));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",0.01000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",3000.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00100000));
+        PyDict_SetItemString(self->params,"sigma_d",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"n_stacking",Py_BuildValue("d",1.00000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CStickyHSStructure.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CStickyHSStructure.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CStickyHSStructure.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,8 +90,8 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"stickiness",Py_BuildValue("d",0.200000));
-        PyDict_SetItemString(self->params,"effect_radius",Py_BuildValue("d",50.000000));
-        PyDict_SetItemString(self->params,"volfraction",Py_BuildValue("d",0.100000));
-        PyDict_SetItemString(self->params,"perturb",Py_BuildValue("d",0.050000));
+        PyDict_SetItemString(self->params,"stickiness",Py_BuildValue("d",0.20000000));
+        PyDict_SetItemString(self->params,"effect_radius",Py_BuildValue("d",50.00000000));
+        PyDict_SetItemString(self->params,"volfraction",Py_BuildValue("d",0.10000000));
+        PyDict_SetItemString(self->params,"perturb",Py_BuildValue("d",0.05000000));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CTriaxialEllipsoidModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CTriaxialEllipsoidModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CTriaxialEllipsoidModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,13 +90,13 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"axis_psi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"semi_axisA",Py_BuildValue("d",35.000000));
-        PyDict_SetItemString(self->params,"semi_axisB",Py_BuildValue("d",100.000000));
-        PyDict_SetItemString(self->params,"semi_axisC",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000005));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"axis_psi",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"semi_axisA",Py_BuildValue("d",35.00000000));
+        PyDict_SetItemString(self->params,"semi_axisB",Py_BuildValue("d",100.00000000));
+        PyDict_SetItemString(self->params,"semi_axisC",Py_BuildValue("d",400.00000000));
+        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.00000530));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/CVesicleModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CVesicleModel.cpp	(revision 308052778784baf067479c9e63914a0293dc7ebd)
+++ sansmodels/src/sans/models/c_models/CVesicleModel.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -90,10 +90,10 @@
         
         // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.000006));
-        PyDict_SetItemString(self->params,"thickness",Py_BuildValue("d",30.000000));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",100.000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"shell_sld",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.00000636));
+        PyDict_SetItemString(self->params,"thickness",Py_BuildValue("d",30.00000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.00000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",100.00000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.00000000));
+        PyDict_SetItemString(self->params,"shell_sld",Py_BuildValue("d",0.00000050));
         // Initialize dispersion / averaging parameter dict
         DispersionVisitor* visitor = new DispersionVisitor();
Index: sansmodels/src/sans/models/c_models/HayterMSA.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/HayterMSA.cpp	(revision f9bf66165686c721c2f0b3a197f7dda1822721ad)
+++ sansmodels/src/sans/models/c_models/HayterMSA.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -36,5 +36,4 @@
 	effect_radius.set_min(0.0);
 	charge      = Parameter(19.0, true);
-	charge.set_min(0.0);
 	volfraction = Parameter(0.0192, true);
 	volfraction.set_min(0.0);
@@ -57,5 +56,5 @@
 	// Add the background after averaging
 	dp[0] = 2.0*effect_radius();
-	dp[1] = charge();
+	dp[1] = fabs(charge());
 	dp[2] = volfraction();
 	dp[3] = temperature();
Index: sansmodels/src/sans/models/c_models/WrapperGenerator.py
===================================================================
--- sansmodels/src/sans/models/c_models/WrapperGenerator.py	(revision 25a608f5e19bd180280be5c566de99df08fb9e8f)
+++ sansmodels/src/sans/models/c_models/WrapperGenerator.py	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -298,5 +298,5 @@
             param_str = "// Initialize parameter dictionary\n"            
             for par in self.params:
-                param_str += "        PyDict_SetItemString(self->params,\"%s\",Py_BuildValue(\"d\",%f));\n" % \
+                param_str += "        PyDict_SetItemString(self->params,\"%s\",Py_BuildValue(\"d\",%10.8f));\n" % \
                     (par, self.params[par])
 
Index: sansmodels/src/sans/models/c_models/c_models.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/c_models.cpp	(revision 5eb91544a3ee42b368ccec3ebcf8e95f75b020a2)
+++ sansmodels/src/sans/models/c_models/c_models.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -15,4 +15,5 @@
 void addCParallelepipedModel(PyObject *module);
 void addCSphereModel(PyObject *module);
+void addCFuzzySphereModel(PyObject *module);
 void addCHardsphereStructure(PyObject *module);
 void addCStickyHSStructure(PyObject *module);
@@ -20,4 +21,5 @@
 void addCHayterMSAStructure(PyObject *module);
 void addCCoreShellModel(PyObject *module);
+void addCCoreFourShellModel(PyObject *module);
 void addCCoreShellCylinderModel(PyObject *module);
 void addCEllipsoidModel(PyObject *module);
@@ -25,4 +27,5 @@
 void addCTriaxialEllipsoidModel(PyObject *module);
 void addCFlexibleCylinderModel(PyObject *module);
+void addCFlexCylEllipXModel(PyObject *module);
 void addCStackedDisksModel(PyObject *module);
 void addCLamellarPSModel(PyObject *module);
@@ -40,4 +43,7 @@
 void addCVesicleModel(PyObject *module);
 void addCBinaryHSModel(PyObject *module);
+//void addCPolymerExclVolModel(PyObject *module);
+void addCFractalModel(PyObject *module);
+void addCPoly_GaussCoil(PyObject *module);
 
 
@@ -213,6 +219,8 @@
 	addCCoreShellCylinderModel(m);
 	addCCoreShellModel(m);
+	addCCoreFourShellModel(m);
 	addCEllipsoidModel(m);
 	addCSphereModel(m);
+	addCFuzzySphereModel(m);
 	addCHardsphereStructure(m);
 	addCStickyHSStructure(m);
@@ -222,4 +230,5 @@
 	addCTriaxialEllipsoidModel(m);
 	addCFlexibleCylinderModel(m);
+	addCFlexCylEllipXModel(m);
 	addCStackedDisksModel(m);
 	addCLamellarPSModel(m);
@@ -242,4 +251,6 @@
 	addCLorentzian(m);
 	addCVesicleModel(m);
-
-}
+	addCPoly_GaussCoil(m);
+	addCFractalModel(m);
+//	addCPolymerExclVolModel(m);
+}
Index: sansmodels/src/sans/models/c_models/corefourshell.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/corefourshell.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/corefourshell.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,239 @@
+/**
+	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
+ *
+ */
+
+#include <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "libSphere.h"
+	#include "corefourshell.h"
+}
+
+CoreFourShellModel :: CoreFourShellModel() {
+	scale      = Parameter(1.0);
+	rad_core     = Parameter(60.0, true);
+	rad_core.set_min(0.0);
+	sld_core   = Parameter(6.4e-6);
+	thick_shell1     = Parameter(10.0, true);
+	thick_shell1.set_min(0.0);
+	sld_shell1   = Parameter(1.0e-6);
+	thick_shell2     = Parameter(10.0, true);
+	thick_shell2.set_min(0.0);
+	sld_shell2   = Parameter(2.0e-6);
+	thick_shell3     = Parameter(10.0, true);
+	thick_shell3.set_min(0.0);
+	sld_shell3   = Parameter(3.0e-6);
+	thick_shell4     = Parameter(10.0, true);
+	thick_shell4.set_min(0.0);
+	sld_shell4   = Parameter(4.0e-6);
+	sld_solv   = Parameter(6.4e-6);
+	background = Parameter(0.001);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * The NIST IGOR library is used for the actual calculation.
+ * @param q: q-value
+ * @return: function value
+ */
+double CoreFourShellModel :: operator()(double q) {
+	double dp[13];
+
+	// Fill parameter array for IGOR library
+	// Add the background after averaging
+	dp[0] = scale();
+	dp[1] = rad_core();
+	dp[2] = sld_core();
+	dp[3] = thick_shell1();
+	dp[4] = sld_shell1();
+	dp[5] = thick_shell2();
+	dp[6] = sld_shell2();
+	dp[7] = thick_shell3();
+	dp[8] = sld_shell3();
+	dp[9] = thick_shell4();
+	dp[10] = sld_shell4();
+	dp[11] = sld_solv();
+	dp[12] = 0.0;
+
+	// Get the dispersion points for the radius
+	vector<WeightPoint> weights_rad;
+	rad_core.get_weights(weights_rad);
+
+	// Get the dispersion points for the thick 1
+	vector<WeightPoint> weights_s1;
+	thick_shell1.get_weights(weights_s1);
+
+	// Get the dispersion points for the thick 2
+	vector<WeightPoint> weights_s2;
+	thick_shell2.get_weights(weights_s2);
+
+	// Get the dispersion points for the thick 3
+	vector<WeightPoint> weights_s3;
+	thick_shell3.get_weights(weights_s3);
+
+	// Get the dispersion points for the thick 4
+	vector<WeightPoint> weights_s4;
+	thick_shell4.get_weights(weights_s4);
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+	double vol = 0.0;
+
+	// Loop over radius weight points
+	for(int i=0; i<weights_rad.size(); i++) {
+		dp[1] = weights_rad[i].value;
+		// Loop over radius weight points
+		for(int j=0; j<weights_s1.size(); j++) {
+			dp[3] = weights_s1[j].value;
+			// Loop over radius weight points
+			for(int k=0; k<weights_s2.size(); k++) {
+				dp[5] = weights_s2[k].value;
+				// Loop over radius weight points
+				for(int l=0; l<weights_s3.size(); l++) {
+					dp[7] = weights_s3[l].value;
+					// Loop over radius weight points
+					for(int m=0; m<weights_s4.size(); m++) {
+						dp[9] = weights_s4[m].value;
+						//Un-normalize FourShell by volume
+						sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
+							* FourShell(dp, q) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3);
+						//Find average volume
+						vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
+							* pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3);
+
+						norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight;
+					}
+				}
+			}
+		}
+	}
+
+	if (vol != 0.0 && norm != 0.0) {
+		//Re-normalize by avg volume
+		sum = sum/(vol/norm);}
+	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 CoreFourShellModel :: operator()(double qx, double qy) {
+	double q = sqrt(qx*qx + qy*qy);
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the sphere
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double CoreFourShellModel :: evaluate_rphi(double q, double phi) {
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double CoreFourShellModel :: calculate_ER() {
+	CoreFourShellParameters dp;
+
+	dp.scale = scale();
+	dp.rad_core = rad_core();
+	dp.sld_core = sld_core();
+	dp.thick_shell1 = thick_shell1();
+	dp.sld_shell1 = sld_shell1();
+	dp.thick_shell2 = thick_shell2();
+	dp.sld_shell2 = sld_shell2();
+	dp.thick_shell3 = thick_shell3();
+	dp.sld_shell3 = sld_shell3();
+	dp.thick_shell4 = thick_shell4();
+	dp.sld_shell4 = sld_shell4();
+	dp.sld_solv = sld_solv();
+	dp.background = 0.0;
+
+	// Get the dispersion points for the radius
+	vector<WeightPoint> weights_rad;
+	rad_core.get_weights(weights_rad);
+
+	// Get the dispersion points for the thick 1
+	vector<WeightPoint> weights_s1;
+	thick_shell1.get_weights(weights_s1);
+
+	// Get the dispersion points for the thick 2
+	vector<WeightPoint> weights_s2;
+	thick_shell2.get_weights(weights_s2);
+
+	// Get the dispersion points for the thick 3
+	vector<WeightPoint> weights_s3;
+	thick_shell3.get_weights(weights_s3);
+
+	// Get the dispersion points for the thick 4
+	vector<WeightPoint> weights_s4;
+	thick_shell4.get_weights(weights_s4);
+
+	double rad_out = 0.0;
+	// 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.rad_core = weights_rad[i].value;
+		// Loop over radius weight points
+		for(int j=0; j<weights_s1.size(); j++) {
+			dp.thick_shell1 = weights_s1[j].value;
+			// Loop over radius weight points
+			for(int k=0; k<weights_s2.size(); k++) {
+				dp.thick_shell2 = weights_s2[k].value;
+				// Loop over radius weight points
+				for(int l=0; l<weights_s3.size(); l++) {
+					dp.thick_shell3 = weights_s3[l].value;
+					// Loop over radius weight points
+					for(int m=0; m<weights_s4.size(); m++) {
+						dp.thick_shell4 = weights_s4[m].value;
+						//Un-normalize FourShell by volume
+						sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
+							* (dp.rad_core+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4);
+						norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight;
+					}
+				}
+			}
+		}
+	}
+	if (norm != 0){
+		//return the averaged value
+		rad_out =  sum/norm;}
+	else{
+		//return normal value
+		rad_out = dp.rad_core+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4;}
+
+	return rad_out;
+}
Index: sansmodels/src/sans/models/c_models/flexcyl_ellipX.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/flexcyl_ellipX.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/flexcyl_ellipX.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,203 @@
+/**
+	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
+ *
+ */
+
+#include <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "libCylinder.h"
+	#include "libStructureFactor.h"
+	#include "flexcyl_ellipX.h"
+}
+
+FlexCylEllipXModel :: FlexCylEllipXModel() {
+	scale      = Parameter(1.0);
+	length     = Parameter(1000.0, true);
+	length.set_min(0.0);
+	kuhn_length = Parameter(100.0, true);
+	kuhn_length.set_min(0.0);
+	radius  = Parameter(20.0, true);
+	radius.set_min(0.0);
+	axis_ratio = Parameter(1.5);
+	axis_ratio.set_min(0.0);
+	sldCyl   = Parameter(1.0e-6);
+	sldSolv   = Parameter(6.3e-6);
+	background = Parameter(0.0001);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * The NIST IGOR library is used for the actual calculation.
+ * @param q: q-value
+ * @return: function value
+ */
+double FlexCylEllipXModel :: operator()(double q) {
+	double dp[8];
+
+	// Fill parameter array for IGOR library
+	// Add the background after averaging
+	dp[0] = scale();
+	dp[1] = length();
+	dp[2] = kuhn_length();
+	dp[3] = radius();
+	dp[4] = axis_ratio();
+	dp[5] = sldCyl();
+	dp[6] = sldSolv();
+	dp[7] = 0.0;
+
+	// Get the dispersion points for the length
+	vector<WeightPoint> weights_len;
+	length.get_weights(weights_len);
+
+	// Get the dispersion points for the kuhn_length
+	vector<WeightPoint> weights_kuhn;
+	kuhn_length.get_weights(weights_kuhn);
+
+	// Get the dispersion points for the radius
+	vector<WeightPoint> weights_rad;
+	radius.get_weights(weights_rad);
+
+	// Get the dispersion points for the axis_ratio
+	vector<WeightPoint> weights_ratio;
+	axis_ratio.get_weights(weights_ratio);
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+	double vol = 0.0;
+
+	// Loop over length weight points
+	for(int i=0; i< (int)weights_len.size(); i++) {
+		dp[1] = weights_len[i].value;
+
+		// Loop over kuhn_length weight points
+		for(int j=0; j< (int)weights_kuhn.size(); j++) {
+			dp[2] = weights_kuhn[j].value;
+
+			// Loop over radius weight points
+			for(int k=0; k< (int)weights_rad.size(); k++) {
+				dp[3] = weights_rad[k].value;
+				// Loop over axis_ratio weight points
+				for(int l=0; l< (int)weights_ratio.size(); l++) {
+					dp[4] = weights_ratio[l].value;
+
+					//Un-normalize by volume
+					sum += weights_len[i].weight * weights_kuhn[j].weight*weights_rad[k].weight
+						* weights_ratio[l].weight * FlexCyl_Ellip(dp, q)
+						* (pow(weights_rad[k].value,2.0) * weights_ratio[l].value * weights_len[i].value);
+					//Find weighted volume
+					vol += weights_rad[k].weight * weights_kuhn[j].weight
+						* weights_len[i].weight * weights_ratio[l].weight
+						*pow(weights_rad[k].value,2.0)* weights_ratio[l].weight*weights_len[i].value;
+					norm += weights_len[i].weight * weights_kuhn[j].weight
+						*weights_rad[k].weight* weights_ratio[l].weight;
+				}
+			}
+		}
+	}
+	if (vol != 0.0 && norm != 0.0) {
+		//Re-normalize by avg volume
+		sum = sum/(vol/norm);}
+
+	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 FlexCylEllipXModel :: operator()(double qx, double qy) {
+	double q = sqrt(qx*qx + qy*qy);
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the triaxial ellipsoid
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double FlexCylEllipXModel :: evaluate_rphi(double q, double phi) {
+	//double qx = q*cos(phi);
+	//double qy = q*sin(phi);
+	return (*this).operator()(q);
+}
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double FlexCylEllipXModel :: calculate_ER() {
+	FlexCyl_EllipXParameters dp;
+
+	dp.radius  = radius();
+	dp.length     = length();
+	dp.axis_ratio  = axis_ratio();
+
+	double rad_out = 0.0;
+	double suf_rad = sqrt(dp.radius*dp.radius*dp.axis_ratio );
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+
+	// Get the dispersion points for the total length
+	vector<WeightPoint> weights_length;
+	length.get_weights(weights_length);
+
+	// Get the dispersion points for minor radius
+	vector<WeightPoint> weights_radius ;
+	radius.get_weights(weights_radius);
+
+	// Get the dispersion points for axis ratio = major_radius/minor_radius
+	vector<WeightPoint> weights_ratio ;
+	axis_ratio.get_weights(weights_ratio);
+
+	// Loop over major shell weight points
+	for(int i=0; i< (int)weights_length.size(); i++) {
+		dp.length = weights_length[i].value;
+		for(int k=0; k< (int)weights_radius.size(); k++) {
+			dp.radius = weights_radius[k].value;
+			// Loop over axis_ratio weight points
+			for(int l=0; l< (int)weights_ratio.size(); l++) {
+				dp.axis_ratio = weights_ratio[l].value;
+				suf_rad = sqrt(dp.radius  * dp.radius  * dp.axis_ratio);
+				//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+				sum +=weights_length[i].weight * weights_radius[k].weight
+						* weights_ratio[l].weight *DiamCyl(dp.length,suf_rad)/2.0;
+				norm += weights_length[i].weight* weights_radius[k].weight* weights_ratio[l].weight;
+			}
+		}
+	}
+	if (norm != 0){
+		//return the averaged value
+		rad_out =  sum/norm;}
+	else{
+		//return normal value
+		//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+		rad_out = DiamCyl(dp.length,suf_rad)/2.0;}
+
+	return rad_out;
+}
Index: sansmodels/src/sans/models/c_models/flexiblecylinder.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/flexiblecylinder.cpp	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/flexiblecylinder.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -42,6 +42,6 @@
 	radius  = Parameter(20.0, true);
 	radius.set_min(0.0);
-	sldCyl   = Parameter(6.3e-6);
-	sldSolv   = Parameter(1.0e-6);
+	sldCyl   = Parameter(1.0e-6);
+	sldSolv   = Parameter(6.3e-6);
 	background = Parameter(0.0001);
 }
@@ -97,10 +97,10 @@
 				sum += weights_len[i].weight
 					* weights_kuhn[j].weight*weights_rad[k].weight * FlexExclVolCyl(dp, q)
-					* pow(weights_rad[k].value,2)*weights_len[i].value;
+					* pow(weights_rad[k].value,2.0)*weights_len[i].value;
 				//Find average volume
 				vol += weights_rad[k].weight
 					* weights_len[i].weight
 					* weights_kuhn[j].weight
-					*pow(weights_rad[k].value,2)*weights_len[i].value;
+					*pow(weights_rad[k].value,2.0)*weights_len[i].value;
 				norm += weights_len[i].weight
 					* weights_kuhn[j].weight*weights_rad[k].weight;
Index: sansmodels/src/sans/models/c_models/fractal.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/fractal.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/fractal.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,120 @@
+/**
+	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
+ *
+ */
+
+#include <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "libTwoPhase.h"
+	#include "fractal.h"
+}
+
+FractalModel :: FractalModel() {
+	scale      = Parameter(1.0);
+    radius = Parameter(5.0, true);
+    radius.set_min(0.0);
+    fractal_dim = Parameter(2.0);
+    cor_length = Parameter(100.0);
+    sldBlock = Parameter(2.0e-6);
+    sldSolv= Parameter(6.35e-6);
+	background = Parameter(0.0);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * The NIST IGOR library is used for the actual calculation.
+ * @param q: q-value
+ * @return: function value
+ */
+double FractalModel :: operator()(double q) {
+	double dp[7];
+
+	// Fill parameter array for IGOR library
+	// Add the background after averaging
+	dp[0] = scale();
+	dp[1] = radius();
+	dp[2] = fractal_dim();
+	dp[3] = cor_length();
+	dp[4] = sldBlock();
+	dp[5] = sldSolv();
+	dp[6] = 0.0;
+
+	// Get the dispersion points for the radius
+	vector<WeightPoint> weights_rad;
+	radius.get_weights(weights_rad);
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+	double vol = 0.0;
+
+	// Loop over radius weight points
+	for(int i=0; i<weights_rad.size(); i++) {
+		dp[1] = weights_rad[i].value;
+
+		//Un-normalize Fractal by volume
+		sum += weights_rad[i].weight
+			* Fractal(dp, fabs(q)) * pow(weights_rad[i].value,3);
+		//Find average volume
+		vol += weights_rad[i].weight
+			* pow(weights_rad[i].value,3);
+
+		norm += weights_rad[i].weight;
+	}
+
+	if (vol != 0.0 && norm != 0.0) {
+		//Re-normalize by avg volume
+		sum = sum/(vol/norm);}
+	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 FractalModel :: operator()(double qx, double qy) {
+	double q = sqrt(qx*qx + qy*qy);
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the FractalModel
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double FractalModel :: evaluate_rphi(double q, double phi) {
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double FractalModel :: calculate_ER() {
+	//NOT implemented yet!!! 'cause None shape Model
+}
Index: sansmodels/src/sans/models/c_models/fuzzysphere.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/fuzzysphere.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/fuzzysphere.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,152 @@
+/**
+	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 <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "fuzzysphere.h"
+}
+
+FuzzySphereModel :: FuzzySphereModel() {
+	scale      = Parameter(0.01);
+	radius     = Parameter(60.0, true);
+	radius.set_min(0.0);
+	fuzziness  = Parameter(10.0);
+	fuzziness.set_min(0.0);
+	sldSph   = Parameter(1.0e-6);
+	sldSolv   = Parameter(3.0e-6);
+	background = Parameter(0.001);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * The NIST IGOR library is used for the actual calculation.
+ * @param q: q-value
+ * @return: function value
+ */
+double FuzzySphereModel :: operator()(double q) {
+	double dp[6];
+
+	// Fill parameter array for IGOR library
+	// Add the background after averaging
+	dp[0] = scale();
+	dp[1] = radius();
+	dp[2] = fuzziness();
+	dp[3] = sldSph();
+	dp[4] = sldSolv();
+	dp[5] = 0.0;
+
+	// Get the dispersion points for the radius
+	vector<WeightPoint> weights_radius;
+	radius.get_weights(weights_radius);
+
+	// Get the dispersion points for the fuzziness
+	vector<WeightPoint> weights_fuzziness;
+	fuzziness.get_weights(weights_fuzziness);
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+	double norm_vol = 0.0;
+	double vol = 0.0;
+
+	// Loop over radius weight points
+	for(int i=0; i<weights_radius.size(); i++) {
+		dp[1] = weights_radius[i].value;
+		// Loop over fuzziness weight points
+		for(int j=0; j<weights_fuzziness.size(); j++) {
+			dp[2] = weights_fuzziness[j].value;
+
+			//Un-normalize SphereForm by volume
+			sum += weights_radius[i].weight * weights_fuzziness[j].weight
+				* fuzzysphere_kernel(dp, q) * pow(weights_radius[i].value,3);
+			//Find average volume : Note the fuzziness has no contribution to the volume
+			vol += weights_radius[i].weight
+				* pow(weights_radius[i].value,3);
+
+			norm += weights_radius[i].weight * weights_fuzziness[j].weight;
+			norm_vol += weights_radius[i].weight;
+		}
+	}
+
+	if (vol != 0.0 && norm_vol != 0.0) {
+		//Re-normalize by avg volume
+		sum = sum/(vol/norm_vol);}
+	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 FuzzySphereModel :: operator()(double qx, double qy) {
+	double q = sqrt(qx*qx + qy*qy);
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the sphere
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double FuzzySphereModel :: evaluate_rphi(double q, double phi) {
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double FuzzySphereModel :: calculate_ER() {
+	FuzzySphereParameters dp;
+	dp.scale = scale();
+	dp.radius = radius();
+	dp.fuzziness = fuzziness();
+	dp.sldSph = sldSph();
+	dp.sldSolv = sldSolv();
+	dp.background = background();
+	double rad_out = 0.0;
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+
+	// Get the dispersion points for the radius
+	// No need to consider the fuzziness.
+	vector<WeightPoint> weights_radius;
+	radius.get_weights(weights_radius);
+	// Loop over radius weight points to average the radius value
+	for(int i=0; i<weights_radius.size(); i++) {
+		sum += weights_radius[i].weight
+			* weights_radius[i].value;
+		norm += weights_radius[i].weight;
+	}
+	if (norm != 0){
+		//return the averaged value
+		rad_out =  sum/norm;}
+	else{
+		//return normal value
+		rad_out = radius();}
+
+	return rad_out;
+}
Index: sansmodels/src/sans/models/c_models/models.hh
===================================================================
--- sansmodels/src/sans/models/c_models/models.hh	(revision f10063e8f3c53d2be93abaea9c0d3d4f3061d444)
+++ sansmodels/src/sans/models/c_models/models.hh	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -22,4 +22,6 @@
 	#include "lamellarPS.h"
 	#include "lamellar.h"
+	#include "fuzzysphere.h"
+	#include "polyexclvol.h"
 }
 
@@ -85,4 +87,24 @@
 	// Constructor
 	SphereModel();
+
+	// Operators to get I(Q)
+	double operator()(double q);
+	double operator()(double qx, double qy);
+	double calculate_ER();
+	double evaluate_rphi(double q, double phi);
+};
+
+class FuzzySphereModel{
+public:
+	// Model parameters
+	Parameter radius;
+	Parameter scale;
+	Parameter fuzziness;
+	Parameter sldSph;
+	Parameter sldSolv;
+	Parameter background;
+
+	// Constructor
+	FuzzySphereModel();
 
 	// Operators to get I(Q)
@@ -218,4 +240,31 @@
 };
 
+class CoreFourShellModel{
+public:
+	// Model parameters
+	Parameter scale;
+	Parameter rad_core;
+	Parameter sld_core;
+	Parameter thick_shell1;
+	Parameter sld_shell1;
+	Parameter thick_shell2;
+	Parameter sld_shell2;
+	Parameter thick_shell3;
+	Parameter sld_shell3;
+	Parameter thick_shell4;
+	Parameter sld_shell4;
+	Parameter sld_solv;
+	Parameter background;
+
+	// Constructor
+	CoreFourShellModel();
+
+	// Operators to get I(Q)
+	double operator()(double q);
+	double operator()(double qx, double qy);
+	double calculate_ER();
+	double evaluate_rphi(double q, double phi);
+};
+
 class CoreShellCylinderModel{
 public:
@@ -320,9 +369,29 @@
 	Parameter sldSolv;
 	Parameter background;
-	Parameter axis_theta;
-	Parameter axis_phi;
 
 	// Constructor
 	FlexibleCylinderModel();
+
+	// Operators to get I(Q)
+	double operator()(double q);
+	double operator()(double qx, double qy);
+	double calculate_ER();
+	double evaluate_rphi(double q, double phi);
+};
+
+class FlexCylEllipXModel{
+public:
+	// Model parameters
+	Parameter scale;
+	Parameter length;
+	Parameter kuhn_length;
+	Parameter radius;
+	Parameter axis_ratio;
+	Parameter sldCyl;
+	Parameter sldSolv;
+	Parameter background;
+
+	// Constructor
+	FlexCylEllipXModel();
 
 	// Operators to get I(Q)
@@ -583,3 +652,61 @@
 };
 
+class Poly_GaussCoil{
+public:
+	// Model parameters
+	Parameter rg;
+	Parameter scale;
+	Parameter poly_m;
+	Parameter background;
+
+	// Constructor
+	Poly_GaussCoil();
+
+	// Operators to get I(Q)
+	double operator()(double q);
+	double operator()(double qx, double qy);
+	double calculate_ER();
+	double evaluate_rphi(double q, double phi);
+};
+
+class FractalModel{
+public:
+	// Model parameters
+	Parameter radius;
+	Parameter scale;
+	Parameter fractal_dim;
+	Parameter cor_length;
+	Parameter sldBlock;
+	Parameter sldSolv;
+	Parameter background;
+
+	// Constructor
+	FractalModel();
+
+	// Operators to get I(Q)
+	double operator()(double q);
+	double operator()(double qx, double qy);
+	double calculate_ER();
+	double evaluate_rphi(double q, double phi);
+};
+
+/***
+class PolymerExclVolModel{
+public:
+	// Model parameters
+	Parameter rg;
+	Parameter scale;
+	Parameter m;
+	Parameter background;
+
+	// Constructor
+	PolymerExclVolModel();
+
+	// Operators to get I(Q)
+	double operator()(double q);
+	double operator()(double qx, double qy);
+	double calculate_ER();
+	double evaluate_rphi(double q, double phi);
+};
+***/
 #endif
Index: sansmodels/src/sans/models/c_models/polygausscoil.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/polygausscoil.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
+++ sansmodels/src/sans/models/c_models/polygausscoil.cpp	(revision c724ccd87499a21402fee0d11d517066f99ddf9b)
@@ -0,0 +1,147 @@
+/**
+	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
+ *
+ */
+
+#include <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "libTwoPhase.h"
+	#include "polygausscoil.h"
+}
+
+Poly_GaussCoil :: Poly_GaussCoil() {
+	scale      = Parameter(1.0);
+	rg     = Parameter(60.0, true);
+	rg.set_min(0.0);
+	poly_m   = Parameter(2.0);
+	background = Parameter(0.001);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * The NIST IGOR library is used for the actual calculation.
+ * @param q: q-value
+ * @return: function value
+ */
+double Poly_GaussCoil :: operator()(double q) {
+	double dp[4];
+
+	// Fill parameter array for IGOR library
+	// Add the background after averaging
+	dp[0] = scale();
+	dp[1] = rg();
+	dp[2] = poly_m();
+	dp[3] = 0.0;
+
+	// Get the dispersion points for the radius
+	//vector<WeightPoint> weights_rad;
+	//rg.get_weights(weights_rad);
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+
+	/**** No rg polydispersion is applied
+	double norm = 0.0;
+	double vol = 0.0;
+
+	// Loop over radius weight points
+	for(int i=0; i<weights_rad.size(); i++) {
+		dp[1] = weights_rad[i].value;
+
+		//Un-normalize SphereForm by volume
+		sum += weights_rad[i].weight
+			* PolyGaussCoil(dp, q) * pow(weights_rad[i].value,3);
+		//Find average volume
+		vol += weights_rad[i].weight
+			* pow(weights_rad[i].value,3);
+
+		norm += weights_rad[i].weight;
+	}
+
+	if (vol != 0.0 && norm != 0.0) {
+		//Re-normalize by avg volume
+		sum = sum/(vol/norm);}
+	return sum/norm + background();
+	****/
+
+	sum = PolyGaussCoil(dp, q);
+	return sum + 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 Poly_GaussCoil :: operator()(double qx, double qy) {
+	double q = sqrt(qx*qx + qy*qy);
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the sphere
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double Poly_GaussCoil :: evaluate_rphi(double q, double phi) {
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double Poly_GaussCoil :: calculate_ER() {
+	PolyGaussCoilParameters dp;
+	dp.scale = scale();
+	dp.rg= rg();
+	dp.poly_m = poly_m();
+	dp.background = background();
+	double rad_out = 0.0;
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+
+	// Get the dispersion points for the radius
+	vector<WeightPoint> weights_rad;
+	rg.get_weights(weights_rad);
+	// Loop over radius weight points to average the radius value
+	for(int i=0; i<weights_rad.size(); i++) {
+		sum += weights_rad[i].weight
+			* weights_rad[i].value;
+		norm += weights_rad[i].weight;
+	}
+	if (norm != 0){
+		//return the averaged value
+		rad_out =  sum/norm;}
+	else{
+		//return normal value
+		rad_out = rg();}
+
+	return rad_out;
+}
