Index: sansmodels/src/sans/models/ParallelepipedModel.py
===================================================================
--- sansmodels/src/sans/models/ParallelepipedModel.py	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
+++ sansmodels/src/sans/models/ParallelepipedModel.py	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+"""
+	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
+"""
+
+""" Provide functionality for a C extension model
+
+	WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ 	         DO NOT MODIFY THIS FILE, MODIFY ..\c_extensions\parallelepiped.h
+ 	         AND RE-RUN THE GENERATOR SCRIPT
+
+"""
+
+from sans.models.BaseComponent import BaseComponent
+from sans_extension.c_models import CParallelepipedModel
+import copy    
+    
+class ParallelepipedModel(CParallelepipedModel, BaseComponent):
+    """ Class that evaluates a ParallelepipedModel model. 
+    	This file was auto-generated from ..\c_extensions\parallelepiped.h.
+    	Refer to that file and the structure it contains
+    	for details of the model.
+    	List of default parameters:
+         scale           = 1.0 
+         short_edgeA     = 35.0 [A]
+         longer_edgeB    = 75.0 [A]
+         longuest_edgeC  = 400.0 [A]
+         contrast        = 5.3e-006 [1/A²]
+         background      = 0.0 [1/cm]
+         parallel_theta  = 1.0 [rad]
+         parallel_phi    = 1.0 [rad]
+
+    """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        BaseComponent.__init__(self)
+        CParallelepipedModel.__init__(self)
+        
+        ## Name of the model
+        self.name = "ParallelepipedModel"
+        ## Model description
+        self.description =""" Calculates the form factor for a rectangular solid with uniform scattering length density.
+		
+		scale:Scale factor
+		short_edgeA: Shortest edge of the parallelepiped [A]
+		longer_edgeB: Longer edge of the parallelepiped [A]
+		longuest_edgeC: Longuest edge of the parallelepiped [A]
+		constrast: particle_sld - solvent_sld
+		background:Incoherent Background [1/cm]"""
+       
+		## Parameter details [units, min, max]
+        self.details = {}
+        self.details['scale'] = ['', None, None]
+        self.details['short_edgeA'] = ['[A]', None, None]
+        self.details['longer_edgeB'] = ['[A]', None, None]
+        self.details['longuest_edgeC'] = ['[A]', None, None]
+        self.details['contrast'] = ['[1/A²]', None, None]
+        self.details['background'] = ['[1/cm]', None, None]
+        self.details['parallel_theta'] = ['[rad]', None, None]
+        self.details['parallel_phi'] = ['[rad]', None, None]
+
+		## fittable parameters
+        self.fixed=['short_edgeA.width', 'longer_edgeB.width', 'longuest_edgeC.width']
+        
+        ## parameters with orientation
+        self.orientation_params =['parallel_phi', 'parallel_theta', 'parallel_phi.width', 'parallel_theta.width']
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        return self._clone(ParallelepipedModel())   
+   
+    def run(self, x = 0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        
+        return CParallelepipedModel.run(self, x)
+   
+    def runXY(self, x = 0.0):
+        """ Evaluate the model in cartesian coordinates
+            @param x: input q, or [qx, qy]
+            @return: scattering function P(q)
+        """
+        
+        return CParallelepipedModel.runXY(self, x)
+        
+    def set_dispersion(self, parameter, dispersion):
+        """
+            Set the dispersion object for a model parameter
+            @param parameter: name of the parameter [string]
+            @dispersion: dispersion object of type DispersionModel
+        """
+        return CParallelepipedModel.set_dispersion(self, parameter, dispersion.cdisp)
+        
+   
+# End of file
Index: sansmodels/src/sans/models/c_extensions/CParallelepipedModel.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/CParallelepipedModel.c	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
+++ sansmodels/src/sans/models/c_extensions/CParallelepipedModel.c	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
@@ -0,0 +1,305 @@
+/** CParallelepipedModel
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY ..\c_extensions\parallelepiped.h
+ *          AND RE-RUN THE GENERATOR SCRIPT
+ *
+ * @author   M.Doucet / UTK
+ */
+ 
+#include <Python.h>
+#include "structmember.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+
+#include "..\c_extensions\parallelepiped.h"
+
+/// Error object for raised exceptions
+static PyObject * CParallelepipedModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Log for unit testing
+    PyObject * log;
+    /// Model parameters
+	ParallelepipedParameters model_pars;
+} CParallelepipedModel;
+
+
+static void
+CParallelepipedModel_dealloc(CParallelepipedModel* self)
+{
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CParallelepipedModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CParallelepipedModel *self;
+    
+    self = (CParallelepipedModel *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CParallelepipedModel_init(CParallelepipedModel *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"longer_edgeB",Py_BuildValue("d",75.000000));
+        PyDict_SetItemString(self->params,"longuest_edgeC",Py_BuildValue("d",400.000000));
+        PyDict_SetItemString(self->params,"parallel_phi",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"parallel_theta",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"short_edgeA",Py_BuildValue("d",35.000000));
+        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000005));
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CParallelepipedModel_members[] = {
+    {"params", T_OBJECT, offsetof(CParallelepipedModel, params), 0,
+     "Parameters"},
+    {"log", T_OBJECT, offsetof(CParallelepipedModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CParallelepipedModel_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 q or [q,phi]
+ * @return: function value
+ */
+static PyObject * run(CParallelepipedModel *self, PyObject *args) {
+	double q_value, phi_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	// Reader parameter dictionary
+    self->model_pars.scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model_pars.longer_edgeB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longer_edgeB") );
+    self->model_pars.longuest_edgeC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longuest_edgeC") );
+    self->model_pars.parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") );
+    self->model_pars.parallel_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_theta") );
+    self->model_pars.background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model_pars.short_edgeA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "short_edgeA") );
+    self->model_pars.contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") );
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CParallelepipedModelError, 
+	    	"CParallelepipedModel.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(CParallelepipedModelError, 
+	    		"CParallelepipedModel.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 = CParallelepipedModel_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CParallelepipedModel_readDouble(PyList_GET_ITEM(pars,1));
+	    // Skip zero
+	    if (q_value==0) {
+	    	return Py_BuildValue("d",0.0);
+	    }
+		return Py_BuildValue("d",_analytical_2D(&(self->model_pars),q_value,phi_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CParallelepipedModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",_analytical_1D(&(self->model_pars),q_value));
+	}	
+}
+
+/**
+ * Function to call to evaluate model in cartesian coordinates
+ * @param args: input q or [qx, qy]]
+ * @return: function value
+ */
+static PyObject * runXY(CParallelepipedModel *self, PyObject *args) {
+	double qx_value, qy_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	// Reader parameter dictionary
+    self->model_pars.scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model_pars.longer_edgeB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longer_edgeB") );
+    self->model_pars.longuest_edgeC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longuest_edgeC") );
+    self->model_pars.parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") );
+    self->model_pars.parallel_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_theta") );
+    self->model_pars.background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model_pars.short_edgeA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "short_edgeA") );
+    self->model_pars.contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") );
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CParallelepipedModelError, 
+	    	"CParallelepipedModel.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(CParallelepipedModelError, 
+	    		"CParallelepipedModel.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 = CParallelepipedModel_readDouble(PyList_GET_ITEM(pars,0));
+	    qy_value = CParallelepipedModel_readDouble(PyList_GET_ITEM(pars,1));
+		return Py_BuildValue("d",_analytical_2DXY(&(self->model_pars),qx_value,qy_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		qx_value = CParallelepipedModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",_analytical_1D(&(self->model_pars),qx_value));
+	}	
+}
+
+static PyObject * reset(CParallelepipedModel *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+
+static PyMethodDef CParallelepipedModel_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+    {"runXY",      (PyCFunction)runXY     , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy"},
+    {"reset",    (PyCFunction)reset   , METH_VARARGS,
+      "Reset pair correlation"},
+    //{"numerical_1D",      (PyCFunction)numerical_1D     , METH_VARARGS,
+    //  "Evaluate the 1D model at a given Q"},
+   {NULL}
+};
+
+static PyTypeObject CParallelepipedModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CParallelepipedModel",             /*tp_name*/
+    sizeof(CParallelepipedModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CParallelepipedModel_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*/
+    "CParallelepipedModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CParallelepipedModel_methods,             /* tp_methods */
+    CParallelepipedModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CParallelepipedModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CParallelepipedModel_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 addCParallelepipedModel(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CParallelepipedModelType) < 0)
+        return;
+
+    Py_INCREF(&CParallelepipedModelType);
+    PyModule_AddObject(module, "CParallelepipedModel", (PyObject *)&CParallelepipedModelType);
+    
+    d = PyModule_GetDict(module);
+    CParallelepipedModelError = PyErr_NewException("CParallelepipedModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CParallelepipedModelError", CParallelepipedModelError);
+}
+
Index: sansmodels/src/sans/models/c_extensions/c_models.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/c_models.c	(revision db1dcef872729a4a27b2bce37df322d93fcf8eda)
+++ sansmodels/src/sans/models/c_extensions/c_models.c	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
@@ -27,4 +27,5 @@
 
 	addCCylinderModel(m);
+	addCParallelepipedModel(m);
 	addCCoreShellCylinderModel(m);
 	addCCoreShellModel(m);
Index: sansmodels/src/sans/models/c_extensions/parallelepiped.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/parallelepiped.c	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
+++ sansmodels/src/sans/models/c_extensions/parallelepiped.c	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
@@ -0,0 +1,123 @@
+/**
+ * Scattering model for a parallelepiped
+ * TODO: Add 2D analysis
+ */
+
+#include "parallelepiped.h"
+#include <math.h>
+#include "libCylinder.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+
+/**
+ * Function to evaluate 1D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+double parallelepiped_analytical_1D(ParallelepipedParameters *pars, double q) {
+	double dp[6];
+	
+	// Fill paramater array
+	dp[0] = pars->scale;
+	dp[1] = pars->short_edgeA;
+	dp[2] = pars->longer_edgeB;
+	dp[3] = pars->longuest_edgeC;
+	dp[4] = pars->contrast;
+	dp[5] = pars->background;
+
+	// Call library function to evaluate model
+	return Parallelepiped(dp, q);	
+}
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+double parallelepiped_analytical_2DXY(ParallelepipedParameters *pars, double qx, double qy) {
+	double q;
+	q = sqrt(qx*qx+qy*qy);
+    return parallelepiped_analytical_2D_scaled(pars, q, qx/q, qy/q);
+} 
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the Parallelepiped
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double parallelepiped_analytical_2D(ParallelepipedParameters *pars, double q, double phi) {
+    return parallelepiped_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
+} 
+        
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the parallelepiped
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y) {
+	double parallel_x, parallel_y, parallel_z;
+	double q_z;
+	double alpha, vol, cos_val;
+	double aa, mu, uu;
+	double answer;
+	
+	
+
+    // parallelepiped orientation
+    parallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi);
+    parallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi);
+    parallel_z = cos(pars->parallel_theta);
+     
+    // q vector
+    q_z = 0;
+        
+    // Compute the angle btw vector q and the
+    // axis of the parallelepiped
+    cos_val = parallel_x*q_x + parallel_y*q_y + parallel_z*q_z;
+    
+    // The following test should always pass
+    if (fabs(cos_val)>1.0) {
+    	printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+     	return 0;
+    }
+    
+    // Note: cos(alpha) = 0 and 1 will get an
+    // undefined value from PPKernel
+	alpha = acos( cos_val );
+
+	aa = pars->short_edgeA/pars->longer_edgeB;
+	mu = 1.0;
+	uu = 1.0;
+	
+	// Call the IGOR library function to get the kernel
+	answer = PPKernel( aa, mu, uu);
+	
+	// Multiply by contrast^2
+	answer *= pars->contrast*pars->contrast;
+	
+	//normalize by cylinder volume
+	//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel
+    vol = pars->short_edgeA* pars->longer_edgeB * pars->longuest_edgeC;
+	answer *= vol;
+	
+	//convert to [cm-1]
+	answer *= 1.0e8;
+	
+	//Scale
+	answer *= pars->scale;
+	
+	// add in the background
+	answer += pars->background;
+	
+	return answer;
+}
+    
+
Index: sansmodels/src/sans/models/c_extensions/parallelepiped.h
===================================================================
--- sansmodels/src/sans/models/c_extensions/parallelepiped.h	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
+++ sansmodels/src/sans/models/c_extensions/parallelepiped.h	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
@@ -0,0 +1,63 @@
+/*
+	TODO: Add 2D model
+*/
+
+#if !defined(parallelepiped_h)
+#define parallelepiped_h
+/** Structure definition for Parallelepiped parameters
+ * [PYTHONCLASS] = ParallelepipedModel
+ * [DISP_PARAMS] = short_edgeA, longer_edgeB, longuest_edgeC,parallel_phi, parallel_theta
+   [DESCRIPTION] = <text> Calculates the form factor for a rectangular solid with uniform scattering length density.
+		
+		scale:Scale factor
+		short_edgeA: Shortest edge of the parallelepiped [A]
+		longer_edgeB: Longer edge of the parallelepiped [A]
+		longuest_edgeC: Longuest edge of the parallelepiped [A]
+		constrast: particle_sld - solvent_sld
+		background:Incoherent Background [1/cm] 
+		</text>
+	[FIXED]= <text>short_edgeA.width; longer_edgeB.width; longuest_edgeC.width</text>
+	[ORIENTATION_PARAMS]= <text>parallel_phi; parallel_theta; parallel_phi.width; parallel_theta.width</text>
+
+
+ **/
+typedef struct {
+    /// Scale factor
+    //  [DEFAULT]=scale=1.0
+    double scale;
+    /// Shortest edge of the parallelepiped [A]
+    //  [DEFAULT]=short_edgeA=35 [A]
+    double short_edgeA;
+	/// Longer edge of the parallelepiped [A]
+    //  [DEFAULT]=longer_edgeB=75 [A]
+    double longer_edgeB;
+	/// Longuest edge of the parallelepiped [A]
+    //  [DEFAULT]=longuest_edgeC=400 [A]
+    double longuest_edgeC;
+    /// Contrast [1/A²]
+    //  [DEFAULT]=contrast=5.3e-6 [1/A²]
+    double contrast;
+	/// Incoherent Background [1/cm] 
+	//  [DEFAULT]=background=0.0 [1/cm]
+	double background;
+    /// Orientation of the parallelepiped axis w/respect incoming beam [rad]
+    //  [DEFAULT]=parallel_theta=1.0 [rad]
+    double parallel_theta;
+    /// Orientation of the parallelepiped in the plane of the detector [rad]
+    //  [DEFAULT]=parallel_phi=1.0 [rad]
+    double parallel_phi;
+
+
+} ParallelepipedParameters;
+
+
+
+/// 1D scattering function
+double parallelepiped_analytical_1D(ParallelepipedParameters *pars, double q);
+
+/// 2D scattering function
+double parallelepiped_analytical_2D(ParallelepipedParameters *pars, double q, double phi);
+double parallelepiped_analytical_2DXY(ParallelepipedParameters *pars, double qx, double qy);
+double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y);
+
+#endif
Index: sansmodels/src/sans/models/c_models/CParallelepipedModel.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/CParallelepipedModel.cpp	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
+++ sansmodels/src/sans/models/c_models/CParallelepipedModel.cpp	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
@@ -0,0 +1,411 @@
+/**
+	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
+ */
+
+/** CParallelepipedModel
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY parallelepiped.h
+ *          AND RE-RUN THE GENERATOR SCRIPT
+ *
+ */
+ 
+extern "C" {
+#include <Python.h>
+#include "structmember.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "parallelepiped.h"
+}
+
+#include "models.hh"
+#include "dispersion_visitor.hh"
+
+/// Error object for raised exceptions
+static PyObject * CParallelepipedModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Dispersion parameters
+    PyObject * dispersion;
+    /// Underlying model object
+    ParallelepipedModel * model;
+    /// Log for unit testing
+    PyObject * log;
+} CParallelepipedModel;
+
+
+static void
+CParallelepipedModel_dealloc(CParallelepipedModel* self)
+{
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CParallelepipedModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CParallelepipedModel *self;
+    
+    self = (CParallelepipedModel *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CParallelepipedModel_init(CParallelepipedModel *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        self->dispersion = PyDict_New();
+        self->model = new ParallelepipedModel();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"longer_edgeB",Py_BuildValue("d",75.000000));
+        PyDict_SetItemString(self->params,"longuest_edgeC",Py_BuildValue("d",400.000000));
+        PyDict_SetItemString(self->params,"parallel_phi",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"parallel_theta",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"short_edgeA",Py_BuildValue("d",35.000000));
+        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000005));
+        // Initialize dispersion / averaging parameter dict
+        DispersionVisitor* visitor = new DispersionVisitor();
+        PyObject * disp_dict;
+        disp_dict = PyDict_New();
+        self->model->short_edgeA.dispersion->accept_as_source(visitor, self->model->short_edgeA.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "short_edgeA", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->longer_edgeB.dispersion->accept_as_source(visitor, self->model->longer_edgeB.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "longer_edgeB", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->longuest_edgeC.dispersion->accept_as_source(visitor, self->model->longuest_edgeC.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "longuest_edgeC", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->parallel_phi.dispersion->accept_as_source(visitor, self->model->parallel_phi.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "parallel_phi", disp_dict);
+        disp_dict = PyDict_New();
+        self->model->parallel_theta.dispersion->accept_as_source(visitor, self->model->parallel_theta.dispersion, disp_dict);
+        PyDict_SetItemString(self->dispersion, "parallel_theta", disp_dict);
+
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CParallelepipedModel_members[] = {
+    {"params", T_OBJECT, offsetof(CParallelepipedModel, params), 0,
+     "Parameters"},
+	{"dispersion", T_OBJECT, offsetof(CParallelepipedModel, dispersion), 0,
+	  "Dispersion parameters"},     
+    {"log", T_OBJECT, offsetof(CParallelepipedModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CParallelepipedModel_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 q or [q,phi]
+ * @return: function value
+ */
+static PyObject * run(CParallelepipedModel *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->longer_edgeB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longer_edgeB") );
+    self->model->longuest_edgeC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longuest_edgeC") );
+    self->model->parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") );
+    self->model->parallel_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_theta") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->short_edgeA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "short_edgeA") );
+    self->model->contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "short_edgeA");
+    self->model->short_edgeA.dispersion->accept_as_destination(visitor, self->model->short_edgeA.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "longer_edgeB");
+    self->model->longer_edgeB.dispersion->accept_as_destination(visitor, self->model->longer_edgeB.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "longuest_edgeC");
+    self->model->longuest_edgeC.dispersion->accept_as_destination(visitor, self->model->longuest_edgeC.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_phi");
+    self->model->parallel_phi.dispersion->accept_as_destination(visitor, self->model->parallel_phi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_theta");
+    self->model->parallel_theta.dispersion->accept_as_destination(visitor, self->model->parallel_theta.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(CParallelepipedModelError, 
+	    	"CParallelepipedModel.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(CParallelepipedModelError, 
+	    		"CParallelepipedModel.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 = CParallelepipedModel_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CParallelepipedModel_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 = CParallelepipedModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(q_value));
+	}	
+}
+
+/**
+ * Function to call to evaluate model in cartesian coordinates
+ * @param args: input q or [qx, qy]]
+ * @return: function value
+ */
+static PyObject * runXY(CParallelepipedModel *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->longer_edgeB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longer_edgeB") );
+    self->model->longuest_edgeC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longuest_edgeC") );
+    self->model->parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") );
+    self->model->parallel_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_theta") );
+    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model->short_edgeA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "short_edgeA") );
+    self->model->contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") );
+    // Read in dispersion parameters
+    PyObject* disp_dict;
+    DispersionVisitor* visitor = new DispersionVisitor();
+    disp_dict = PyDict_GetItemString(self->dispersion, "short_edgeA");
+    self->model->short_edgeA.dispersion->accept_as_destination(visitor, self->model->short_edgeA.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "longer_edgeB");
+    self->model->longer_edgeB.dispersion->accept_as_destination(visitor, self->model->longer_edgeB.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "longuest_edgeC");
+    self->model->longuest_edgeC.dispersion->accept_as_destination(visitor, self->model->longuest_edgeC.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_phi");
+    self->model->parallel_phi.dispersion->accept_as_destination(visitor, self->model->parallel_phi.dispersion, disp_dict);
+    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_theta");
+    self->model->parallel_theta.dispersion->accept_as_destination(visitor, self->model->parallel_theta.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(CParallelepipedModelError, 
+	    	"CParallelepipedModel.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(CParallelepipedModelError, 
+	    		"CParallelepipedModel.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 = CParallelepipedModel_readDouble(PyList_GET_ITEM(pars,0));
+	    qy_value = CParallelepipedModel_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 = CParallelepipedModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",(*(self->model))(qx_value));
+	}	
+}
+
+static PyObject * reset(CParallelepipedModel *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+static PyObject * set_dispersion(CParallelepipedModel *self, PyObject *args) {
+	PyObject * disp;
+	const char * par_name;
+
+	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
+	    PyErr_SetString(CParallelepipedModelError,
+	    	"CParallelepipedModel.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, "short_edgeA")) {
+        self->model->short_edgeA.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "longer_edgeB")) {
+        self->model->longer_edgeB.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "longuest_edgeC")) {
+        self->model->longuest_edgeC.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "parallel_phi")) {
+        self->model->parallel_phi.dispersion = dispersion;
+    } else    if (!strcmp(par_name, "parallel_theta")) {
+        self->model->parallel_theta.dispersion = dispersion;
+    } else {
+	    PyErr_SetString(CParallelepipedModelError,
+	    	"CParallelepipedModel.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 CParallelepipedModel_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q or Q, phi"},
+    {"runXY",      (PyCFunction)runXY     , METH_VARARGS,
+      "Evaluate the model at a given Q or Qx, Qy"},
+    {"reset",    (PyCFunction)reset   , METH_VARARGS,
+      "Reset pair correlation"},
+    {"set_dispersion",      (PyCFunction)set_dispersion     , METH_VARARGS,
+      "Set the dispersion model for a given parameter"},
+   {NULL}
+};
+
+static PyTypeObject CParallelepipedModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CParallelepipedModel",             /*tp_name*/
+    sizeof(CParallelepipedModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CParallelepipedModel_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*/
+    "CParallelepipedModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CParallelepipedModel_methods,             /* tp_methods */
+    CParallelepipedModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CParallelepipedModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CParallelepipedModel_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 addCParallelepipedModel(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CParallelepipedModelType) < 0)
+        return;
+
+    Py_INCREF(&CParallelepipedModelType);
+    PyModule_AddObject(module, "CParallelepipedModel", (PyObject *)&CParallelepipedModelType);
+    
+    d = PyModule_GetDict(module);
+    CParallelepipedModelError = PyErr_NewException("CParallelepipedModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CParallelepipedModelError", CParallelepipedModelError);
+}
+
Index: sansmodels/src/sans/models/c_models/c_models.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/c_models.cpp	(revision 4d09267b94c669d6a8811368c69c8d130ef27ffb)
+++ sansmodels/src/sans/models/c_models/c_models.cpp	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
@@ -9,4 +9,5 @@
 
 void addCCylinderModel(PyObject *module);
+void addCParallelepipedModel(PyObject *module);
 void addCSphereModel(PyObject *module);
 void addCHardsphereStructure(PyObject *module);
@@ -148,4 +149,5 @@
 
 	addCCylinderModel(m);
+	addCParallelepipedModel(m);
 	addCCoreShellCylinderModel(m);
 	addCCoreShellModel(m);
Index: sansmodels/src/sans/models/c_models/models.hh
===================================================================
--- sansmodels/src/sans/models/c_models/models.hh	(revision a0782f15b80e58c94f6092aa540cf0c010047197)
+++ sansmodels/src/sans/models/c_models/models.hh	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
@@ -19,4 +19,5 @@
 extern "C" {
 	#include "cylinder.h"
+	#include "parallelepiped.h"
 }
 
@@ -42,4 +43,27 @@
 	double evaluate_rphi(double q, double phi);
 };
+
+class ParallelepipedModel{
+public:
+	// TODO: add 2D 
+	// Model parameters
+	Parameter scale;
+	Parameter short_edgeA;
+	Parameter longer_edgeB;
+	Parameter longuest_edgeC;
+	Parameter contrast;
+	Parameter background;
+	Parameter parallel_theta;
+	Parameter parallel_phi;
+
+	// Constructor
+	ParallelepipedModel();
+
+	// Operators to get I(Q)
+	double operator()(double q);
+	double operator()(double qx, double qy);
+	double evaluate_rphi(double q, double phi);
+};
+
 
 class SphereModel{
Index: sansmodels/src/sans/models/c_models/parallelepiped.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/parallelepiped.cpp	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
+++ sansmodels/src/sans/models/c_models/parallelepiped.cpp	(revision 8a487138bf42fe6104144ee6073b5989aa12bcd6)
@@ -0,0 +1,214 @@
+/**
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	If you use DANSE applications to do scientific research that leads to
+	publication, we ask that you acknowledge the use of the software with the
+	following sentence:
+
+	"This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+	copyright 2008, University of Tennessee
+ */
+
+/**
+ * Scattering model classes
+ * The classes use the IGOR library found in
+ *   sansmodels/src/libigor
+ *
+ *	TODO: refactor so that we pull in the old sansmodels.c_extensions
+ *	TODO: add 2D function
+ */
+
+#include <math.h>
+#include "models.hh"
+#include "parameters.hh"
+#include <stdio.h>
+using namespace std;
+
+extern "C" {
+	#include "libCylinder.h"
+	#include "parallelepiped.h"
+}
+
+ParallelepipedModel :: ParallelepipedModel() {
+	scale      = Parameter(1.0);
+	short_edgeA     = Parameter(35.0, true);
+	short_edgeA.set_max(1.0);
+	longer_edgeB     = Parameter(75.0, true);
+	longer_edgeB.set_min(1.0);
+	longuest_edgeC     = Parameter(400.0, true);
+	longuest_edgeC.set_min(1.0);
+	contrast   = Parameter(53.e-7);
+	background = Parameter(0.0);
+	parallel_theta  = Parameter(0.0, true);
+	parallel_phi    = Parameter(0.0, true);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * The NIST IGOR library is used for the actual calculation.
+ * @param q: q-value
+ * @return: function value
+ */
+double ParallelepipedModel :: operator()(double q) {
+	double dp[5];
+
+	// Fill parameter array for IGOR library
+	// Add the background after averaging
+	dp[0] = scale();
+	dp[1] = short_edgeA();
+	dp[2] = longer_edgeB();
+	dp[3] = longuest_edgeC();
+	dp[4] = contrast();
+	//dp[5] = background();
+	dp[5] = 0.0;
+
+	// Get the dispersion points for the short_edgeA
+	vector<WeightPoint> weights_short_edgeA;
+	short_edgeA.get_weights(weights_short_edgeA);
+	
+	// Get the dispersion points for the longer_edgeB
+	vector<WeightPoint> weights_longer_edgeB;
+	longer_edgeB.get_weights(weights_longer_edgeB);
+
+	// Get the dispersion points for the longuest_edgeC
+	vector<WeightPoint> weights_longuest_edgeC;
+	longuest_edgeC.get_weights(weights_longuest_edgeC);
+
+
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+	
+	// Loop over short_edgeA weight points
+	for(int i=0; i< (int)weights_short_edgeA.size(); i++) {
+		dp[1] = weights_short_edgeA[i].value;
+
+		// Loop over longer_edgeB weight points
+		for(int j=0; j< (int)weights_longer_edgeB.size(); j++) {
+			dp[2] = weights_longer_edgeB[i].value;
+
+			// Loop over longuest_edgeC weight points
+			for(int k=0; k< (int)weights_longuest_edgeC.size(); k++) {
+				dp[3] = weights_longuest_edgeC[j].value;
+
+				sum += weights_short_edgeA[i].weight * weights_longer_edgeB[j].weight 
+					* weights_longuest_edgeC[k].weight * Parallelepiped(dp, q);
+
+				norm += weights_short_edgeA[i].weight
+					 * weights_longer_edgeB[j].weight * weights_longuest_edgeC[k].weight;
+			}
+		}
+	}
+	return sum/norm + background();
+}
+/**
+ * Function to evaluate 2D scattering function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double ParallelepipedModel :: operator()(double qx, double qy) {
+	ParallelepipedParameters dp;
+	// Fill parameter array
+	dp.scale      = scale();
+	dp.short_edgeA   = short_edgeA();
+	dp.longer_edgeB   = longer_edgeB();
+	dp.longuest_edgeC  = longuest_edgeC();
+	dp.contrast   = contrast();
+	dp.background = 0.0;
+	//dp.background = background();
+	dp.parallel_theta  = parallel_theta();
+	dp.parallel_phi    = parallel_phi();
+	
+
+	// Get the dispersion points for the short_edgeA
+	vector<WeightPoint> weights_short_edgeA;
+	short_edgeA.get_weights(weights_short_edgeA);
+
+	// Get the dispersion points for the longer_edgeB
+	vector<WeightPoint> weights_longer_edgeB;
+	longer_edgeB.get_weights(weights_longer_edgeB);
+
+	// Get angular averaging for the longuest_edgeC
+	vector<WeightPoint> weights_longuest_edgeC;
+	longuest_edgeC.get_weights(weights_longuest_edgeC);
+
+	// Get angular averaging for theta
+	vector<WeightPoint> weights_parallel_theta;
+	parallel_theta.get_weights(weights_parallel_theta);
+
+	// Get angular averaging for phi
+	vector<WeightPoint> weights_parallel_phi;
+	parallel_phi.get_weights(weights_parallel_phi);
+
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+
+	// Loop over radius weight points
+	for(int i=0; i< (int)weights_short_edgeA.size(); i++) {
+		dp.short_edgeA = weights_short_edgeA[i].value;
+
+		// Loop over longer_edgeB weight points
+		for(int j=0; j< (int)weights_longer_edgeB.size(); j++) {
+			dp.longer_edgeB = weights_longer_edgeB[j].value;
+
+			// Average over longuest_edgeC distribution
+			for(int k=0; k< (int)weights_longuest_edgeC.size(); k++) {
+				dp.longuest_edgeC = weights_longuest_edgeC[k].value;
+
+				// Average over theta distribution
+				for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
+				dp.parallel_theta = weights_parallel_theta[l].value;
+
+					// Average over phi distribution
+					for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
+						dp.parallel_phi = weights_parallel_phi[m].value;
+
+						double _ptvalue = weights_short_edgeA[i].weight
+							* weights_longer_edgeB[j].weight
+							* weights_longuest_edgeC[k].weight
+							* weights_parallel_theta[l].weight
+							* weights_parallel_phi[m].weight
+							* parallelepiped_analytical_2DXY(&dp, qx, qy);
+						if (weights_parallel_theta.size()>1) {
+							_ptvalue *= sin(weights_parallel_theta[l].value);
+						}
+						sum += _ptvalue;
+
+						norm += weights_short_edgeA[i].weight
+							* weights_longer_edgeB[j].weight
+							* weights_longuest_edgeC[k].weight
+							* weights_parallel_theta[l].weight
+							* weights_parallel_phi[m].weight;
+					}
+
+				}
+			}
+		}
+	}
+	// Averaging in theta needs an extra normalization
+	// factor to account for the sin(theta) term in the
+	// integration (see documentation).
+	if (weights_parallel_theta.size()>1) norm = norm / asin(1.0);
+	return sum/norm + background();
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double ParallelepipedModel :: evaluate_rphi(double q, double phi) {
+	double qx = q*cos(phi);
+	double qy = q*sin(phi);
+	return (*this).operator()(qx, qy);
+}
