Index: sansmodels/src/sans/models/prototypes/DispCylinderModel.py
===================================================================
--- sansmodels/src/sans/models/prototypes/DispCylinderModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/DispCylinderModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+""" Provide functionality for a C extension model
+
+	WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ 	         DO NOT MODIFY THIS FILE, MODIFY disperse_cylinder.h
+ 	         AND RE-RUN THE GENERATOR SCRIPT
+
+    @author: Mathieu Doucet / UTK
+    @contact: mathieu.doucet@nist.gov
+"""
+
+from sans.models.BaseComponent import BaseComponent
+from sans_extension.prototypes.c_models import CDispCylinderModel
+import copy    
+    
+class DispCylinderModel(CDispCylinderModel, BaseComponent):
+    """ Class that evaluates a DispCylinderModel model. 
+    	This file was auto-generated from disperse_cylinder.h.
+    	Refer to that file and the structure it contains
+    	for details of the model.
+    	List of default parameters:
+         scale           = 1.0 
+         radius          = 20.0 A
+         length          = 400.0 A
+         contrast        = 3e-006 A-2
+         background      = 0.0 cm-1
+         cyl_theta       = 1.57 rad
+         cyl_phi         = 0.0 rad
+         cyl_phi         = 0.0 rad
+         sigma_theta     = 0.0 rad
+         sigma_phi       = 0.0 rad
+         sigma_radius    = 0.0 rad
+         sigma_radius    = 0.0 rad
+         n_pts           = 25.0 
+
+    """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        BaseComponent.__init__(self)
+        CDispCylinderModel.__init__(self)
+        
+        ## Name of the model
+        self.name = "DispCylinderModel"
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        obj = DispCylinderModel()
+        obj.params = copy.deepcopy(self.params)
+        return obj   
+   
+    def run(self, x = 0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        
+        return CDispCylinderModel.run(self, x)
+   
+# End of file
Index: sansmodels/src/sans/models/prototypes/DisperseCylinderModel.py
===================================================================
--- sansmodels/src/sans/models/prototypes/DisperseCylinderModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/DisperseCylinderModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+""" Provide functionality for a C extension model
+
+
+    @author: Mathieu Doucet / UTK
+    @contact: mathieu.doucet@nist.gov
+"""
+
+from sans.models.CylinderModel import CylinderModel
+from sans_extension.c_models import Disperser
+from sans.models.DisperseModel import DisperseModel
+import copy    
+    
+class DisperseCylinderModel(CylinderModel):
+    """ Class that evaluates a SmearCylinderModel model. 
+    """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        CylinderModel.__init__(self)
+        self.params['sigma_phi'] = 0.0
+        self.params['sigma_theta'] = 0.0
+        self.params['n_pts'] = 20
+        
+        
+        ## Name of the model
+        self.name = "DisperseCylinderModel"
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        obj = DisperseCylinderModel()
+        obj.params = copy.deepcopy(self.params)
+        return obj   
+   
+    def run(self, x = 0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        model = CylinderModel()
+        model.params = copy.deepcopy(self.params)
+        #d = Disperser(model, ["cyl_phi", "cyl_theta"], 
+		#				  [self.params['sigma_phi'], self.params['sigma_theta']])
+        d = DisperseModel(model, ["cyl_phi", "cyl_theta"], 
+					  [self.params['sigma_phi'], self.params['sigma_theta']])
+        d.params['n_pts'] = self.params['n_pts']
+        return d.run(x)
+   
+# End of file
Index: sansmodels/src/sans/models/prototypes/DisperseEllipsoidModel.py
===================================================================
--- sansmodels/src/sans/models/prototypes/DisperseEllipsoidModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/DisperseEllipsoidModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,50 @@
+#!/usr/bin/env python
+""" Provide functionality for a C extension model
+
+
+    @author: Mathieu Doucet / UTK
+    @contact: mathieu.doucet@nist.gov
+"""
+
+from sans.models.EllipsoidModel import EllipsoidModel
+from sans.models.DisperseModel import DisperseModel
+import copy    
+    
+class DisperseEllipsoidModel(EllipsoidModel):
+    """ Class that evaluates a DisperseEllipsoidModel model. 
+    """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        EllipsoidModel.__init__(self)
+        self.params['sigma_phi'] = 0.0
+        self.params['sigma_theta'] = 0.0
+        self.params['n_pts'] = 20
+        
+        
+        ## Name of the model
+        self.name = "DisperseEllipsoidModel"
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        obj = DisperseEllipsoidModel()
+        obj.params = copy.deepcopy(self.params)
+        return obj   
+   
+    def run(self, x = 0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        model = EllipsoidModel()
+        model.params = copy.deepcopy(self.params)
+        #d = Disperser(model, ["cyl_phi", "cyl_theta"], 
+		#				  [self.params['sigma_phi'], self.params['sigma_theta']])
+        d = DisperseModel(model, ["axis_phi", "axis_theta"], 
+					  [self.params['sigma_phi'], self.params['sigma_theta']])
+        d.params['n_pts'] = self.params['n_pts']
+        return d.run(x)
+   
+# End of file
Index: sansmodels/src/sans/models/prototypes/SimCylinder.py
===================================================================
--- sansmodels/src/sans/models/prototypes/SimCylinder.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/SimCylinder.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+""" Provide functionality for a C extension model
+
+	WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ 	         DO NOT MODIFY THIS FILE, MODIFY simcylinder.h
+ 	         AND RE-RUN THE GENERATOR SCRIPT
+
+    @author: Mathieu Doucet / UTK
+    @contact: mathieu.doucet@nist.gov
+"""
+
+from sans.models.BaseComponent import BaseComponent
+from sans_extension.prototypes.c_models import CSimCylinder
+import copy    
+    
+class SimCylinder(CSimCylinder, BaseComponent):
+    """ Class that evaluates a SimCylinder model. 
+    	This file was auto-generated from simcylinder.h.
+    	Refer to that file and the structure it contains
+    	for details of the model.
+    	List of default parameters:
+         scale           = 1.0 
+         radius          = 20.0 
+         length          = 400.0 
+         theta           = 0.0 
+         phi             = 0.0 
+         qmax            = 0.1 
+
+    """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        BaseComponent.__init__(self)
+        CSimCylinder.__init__(self)
+        
+        ## Name of the model
+        self.name = "SimCylinder"
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        obj = SimCylinder()
+        obj.params = copy.deepcopy(self.params)
+        return obj   
+   
+    def run(self, x = 0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        
+        return CSimCylinder.run(self, x)
+# main
+if __name__ == '__main__':
+    app = SimCylinder()
+    print "%s: f(1) = %g" % (app.name, app.run(1))
+    print "%s: f(1,1.3) = %g" % (app.name, app.run([1, 1.3]))
+   
+# End of file
Index: sansmodels/src/sans/models/prototypes/SimCylinderF.py
===================================================================
--- sansmodels/src/sans/models/prototypes/SimCylinderF.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/SimCylinderF.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,54 @@
+#!/usr/bin/env python
+""" Provide functionality for a C extension model
+
+	WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ 	         DO NOT MODIFY THIS FILE, MODIFY ..\prototypes\src\simcylinder_fast.h
+ 	         AND RE-RUN THE GENERATOR SCRIPT
+
+    @author: Mathieu Doucet / UTK
+    @contact: mathieu.doucet@nist.gov
+"""
+
+from sans.models.BaseComponent import BaseComponent
+from sans_extension.prototypes.c_models import CSimCylinderF
+import copy    
+    
+class SimCylinderF(CSimCylinderF, BaseComponent):
+    """ Class that evaluates a SimCylinderF model. 
+    	This file was auto-generated from ..\prototypes\src\simcylinder_fast.h.
+    	Refer to that file and the structure it contains
+    	for details of the model.
+    	List of default parameters:
+         scale           = 1.0 
+         radius          = 20.0 
+         length          = 400.0 
+         theta           = 0.0 
+         phi             = 0.0 
+
+    """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        BaseComponent.__init__(self)
+        CSimCylinderF.__init__(self)
+        
+        ## Name of the model
+        self.name = "SimCylinderF"
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        obj = SimCylinderF()
+        obj.params = copy.deepcopy(self.params)
+        return obj   
+   
+    def run(self, x = 0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        
+        return CSimCylinderF.run(self, x)
+   
+# End of file
Index: sansmodels/src/sans/models/prototypes/SimSphereF.py
===================================================================
--- sansmodels/src/sans/models/prototypes/SimSphereF.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/SimSphereF.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+""" Provide functionality for a C extension model
+
+	WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ 	         DO NOT MODIFY THIS FILE, MODIFY ../prototypes/src/sphere_fast.h
+ 	         AND RE-RUN THE GENERATOR SCRIPT
+
+    @author: Mathieu Doucet / UTK
+    @contact: mathieu.doucet@nist.gov
+"""
+
+from sans.models.BaseComponent import BaseComponent
+from sans_extension.prototypes.c_models import CSimSphereF
+import copy    
+    
+class SimSphereF(CSimSphereF, BaseComponent):
+    """ Class that evaluates a SimSphereF model. 
+    	This file was auto-generated from ../prototypes/src/sphere_fast.h.
+    	Refer to that file and the structure it contains
+    	for details of the model.
+    	List of default parameters:
+         scale           = 1.0 
+         radius          = 20.0 
+
+    """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        BaseComponent.__init__(self)
+        CSimSphereF.__init__(self)
+        
+        ## Name of the model
+        self.name = "SimSphereF"
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        obj = SimSphereF()
+        obj.params = copy.deepcopy(self.params)
+        return obj   
+   
+    def run(self, x = 0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        
+        return CSimSphereF.run(self, x)
+   
+# End of file
Index: sansmodels/src/sans/models/prototypes/SmearCylinderModel.py
===================================================================
--- sansmodels/src/sans/models/prototypes/SmearCylinderModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/SmearCylinderModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+""" Provide functionality for a C extension model
+
+	WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ 	         DO NOT MODIFY THIS FILE, MODIFY smeared_cylinder.h
+ 	         AND RE-RUN THE GENERATOR SCRIPT
+
+    @author: Mathieu Doucet / UTK
+    @contact: mathieu.doucet@nist.gov
+"""
+
+from sans.models.BaseComponent import BaseComponent
+from sans_extension.prototypes.c_models import CSmearCylinderModel
+import copy    
+    
+class SmearCylinderModel(CSmearCylinderModel, BaseComponent):
+    """ Class that evaluates a SmearCylinderModel model. 
+    	This file was auto-generated from smeared_cylinder.h.
+    	Refer to that file and the structure it contains
+    	for details of the model.
+    	List of default parameters:
+         scale           = 1.0 
+         radius          = 20.0 A
+         length          = 400.0 A
+         contrast        = 3e-006 A-2
+         background      = 0.0 cm-1
+         cyl_theta       = 1.0 rad
+         cyl_phi         = 1.0 rad
+         cyl_phi         = 0.0 rad
+         sigma_theta     = 0.0 rad
+         sigma_phi       = 0.0 rad
+         sigma_radius    = 0.0 rad
+
+    """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        BaseComponent.__init__(self)
+        CSmearCylinderModel.__init__(self)
+        
+        ## Name of the model
+        self.name = "SmearCylinderModel"
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        obj = SmearCylinderModel()
+        obj.params = copy.deepcopy(self.params)
+        return obj   
+   
+    def run(self, x = 0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        
+        return CSmearCylinderModel.run(self, x)
+   
+# End of file
Index: sansmodels/src/sans/models/prototypes/TestSphere2.py
===================================================================
--- sansmodels/src/sans/models/prototypes/TestSphere2.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/TestSphere2.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+""" Provide functionality for a C extension model
+
+	WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ 	         DO NOT MODIFY THIS FILE, MODIFY testsphere.h
+ 	         AND RE-RUN THE GENERATOR SCRIPT
+
+    @author: Mathieu Doucet / UTK
+    @contact: mathieu.doucet@nist.gov
+"""
+
+from sans.models.BaseComponent import BaseComponent
+from sans_extension.prototypes.c_models import CTestSphere2
+import copy    
+    
+class TestSphere2(CTestSphere2, BaseComponent):
+    """ Scattering model based on a cylinder """
+        
+    def __init__(self):
+        """ Initialization """
+        
+        # Initialize BaseComponent first, then sphere
+        BaseComponent.__init__(self)
+        CTestSphere2.__init__(self)
+        
+        ## Name of the model
+        self.name = "TestSphere2"
+   
+    def clone(self):
+        """ Return a identical copy of self """
+        obj = TestSphere2()
+        obj.params = copy.deepcopy(self.params)
+        return obj   
+    
+    def run(self, x=0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        
+        return CTestSphere2.run(self, x)
+   
+# End of file
Index: sansmodels/src/sans/models/prototypes/WeightModel.py
===================================================================
--- sansmodels/src/sans/models/prototypes/WeightModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/WeightModel.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,119 @@
+#!/usr/bin/env python
+""" 
+    Wrapper for the Disperser class extension
+
+    @author: Mathieu Doucet / UTK
+    @contact: mathieu.doucet@nist.gov
+"""
+
+import math
+import os, string, sys
+from sans.models.BaseComponent import BaseComponent
+from sans.models.CylinderModel import CylinderModel
+    
+class WeightModel(CylinderModel):
+    """ 
+    """
+        
+    def __init__(self):
+        """ Initialization 
+            @param model: Model to disperse [BaseComponent]
+            @param param: Parameter to disperse [string]
+        """
+        
+        # Initialize BaseComponent first, then sphere
+        CylinderModel.__init__(self)
+        self.params['datafile'] = '' 
+        
+        self.d = None
+        ## Name of the model
+        self.name = 'WeightModel'
+   
+    def run(self, x = 0.0):
+        """ Evaluate the model
+            @param x: input q, or [q,phi]
+            @return: scattering function P(q)
+        """
+        if len(self.d.x)==0:
+            return 0.0
+        try:
+            total = 0
+            norma = 0
+            for i in range(len(self.d.x)):
+                self.params['cyl_phi'] = self.d.x[i]*math.pi/180.0
+                total += CylinderModel.run(self, x) * self.d.y[i]
+                
+                
+                self.params['cyl_phi'] = -self.d.x[i]*math.pi/180.0
+                total += CylinderModel.run(self, x) * self.d.y[i]
+                
+                norma += 2.0*self.d.y[i]
+            return total/norma
+        except:
+            print sys.exc_value
+            
+        return 0.0
+           
+    def readData(self):
+                self.d = DataReader(self.params['datafile'])
+                self.d.read()
+
+   
+class DataReader:
+    """ Simple data reader for Igor data files """
+    
+    def __init__(self, filename):
+        """ Init
+            @param filename: Name of Igor data file to read
+        """
+        self.file = filename
+        self.x = []
+        self.y = []
+        
+    def read(self):
+        """ Read file """
+        # Check if the file is there
+        if not os.path.isfile(self.file):
+            raise ValueError, \
+            "Specified file %s is not a regular file" % self.file
+        
+        # Read file
+        f = open(self.file,'r')
+        buf = f.read()
+        
+        # Get content
+        dataStarted = False
+        
+        
+        lines = string.split(buf,'\n')
+        itot = 0
+        self.x = []
+        self.y = []
+        for line in lines:
+            if line.count("Angle")>0:
+                dataStarted = True
+                continue
+                
+            if dataStarted == True:
+                
+                try:
+                    toks = line.split()
+                except:
+                    print sys.exc_value
+                
+                self.x.append(float(toks[0]))
+                
+                self.y.append(float(toks[1]))
+                
+                
+                itot += 1
+                    
+        print "Read %g points from file %s" % (len(self.image), self.file)
+                
+       
+if __name__ == '__main__':
+    app = WeightModel()
+    cyl = CylinderModel()
+    print "%s: f(1) = %g" % (app.name, app.run(1))
+   
+# End of file
Index: sansmodels/src/sans/models/prototypes/__init__.py
===================================================================
--- sansmodels/src/sans/models/prototypes/__init__.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/__init__.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,8 @@
+"""
+    Extension for testing of SANS models
+    
+    This module is experimental and serves as a repository for prototypes
+    and tests for the modeling functionality for SANS.
+    
+    Use this module at your own risks.
+"""
Index: sansmodels/src/sans/models/prototypes/setup.py
===================================================================
--- sansmodels/src/sans/models/prototypes/setup.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/setup.py	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,61 @@
+"""
+ Installation script for SANS models
+
+  - To compile and install:
+      python setup.py install
+  - To create distribution:
+      python setup.py bdist_wininst
+
+"""
+
+from distutils.core import setup, Extension
+
+
+# Build the module name
+srcdir  = "src"
+igordir = "../../../libigor"
+extdir  = "../c_extensions"
+
+print "Installing SANS models"
+
+setup(
+    name="models",
+    version = "0.1",
+    description = "Python module for various tests of SANS scattering models",
+    author = "Mathieu Doucet",
+    author_email = "doucet@nist.gov",
+    url = "http://danse.us/trac/sans",
+    
+    # Place this module under the sans package
+    #ext_package = "sans",
+    
+    # Use the pure python modules
+    package_dir = {"sans_extension.prototypes":"."},
+    
+    packages = ["sans_extension.prototypes"],
+    
+    ext_modules = [ Extension("sans_extension.prototypes.c_models",
+     sources = [
+        srcdir+"/c_models.c",
+        "../../../igor_wrapper/c_disperser.c",
+        srcdir+"/disperse_cylinder.c",
+        srcdir+"/CDispCylinderModel.c",
+        srcdir+"/CSimCylinderF.c",
+        srcdir+"/simcylinder_fast.c",
+        srcdir+"/CSimSphereF.c",
+        srcdir+"/sphere_fast.c",
+        srcdir+"/CTestSphere2.c",
+        srcdir+"/testsphere.c",
+        srcdir+"/CSimCylinder.c",
+        srcdir+"/simcylinder.c",
+        srcdir+"/modelCalculations.c",
+        srcdir+"/CSmearCylinderModel.c",
+        srcdir+"/smeared_cylinder.c",
+        igordir+"/libCylinder.c",
+        igordir+"/libSphere.c",
+        extdir+"/cylinder.c",
+        srcdir+"/CCanvas.c",
+        srcdir+"/canvas.c"
+            ],
+         include_dirs=[igordir, extdir])])
+        
Index: sansmodels/src/sans/models/prototypes/src/CCanvas.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/CCanvas.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/CCanvas.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,223 @@
+
+#include <Python.h>
+#include "structmember.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+
+#include "canvas.h"
+
+/// Error object for raised exceptions
+static PyObject * CCanvasError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Model parameters
+	CanvasParams canvas_pars;
+} CCanvas;
+
+
+static void
+CCanvas_dealloc(CCanvas* self)
+{
+    canvas_dealloc(&(self->canvas_pars));
+     
+    self->ob_type->tp_free((PyObject*)self);
+
+}
+
+static PyObject *
+CCanvas_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CCanvas *self;
+    
+    self = (CCanvas *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CCanvas_init(CCanvas *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+		srand(1001);	
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+  		canvas_init(&(self->canvas_pars));
+    }
+    return 0;
+}
+
+static PyMemberDef CCanvas_members[] = {
+    {"params", T_OBJECT, offsetof(CCanvas, params), 0,
+     "Parameters"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CCanvas_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(CCanvas *self, PyObject *args) {
+	double q_value, phi_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CCanvasError, 
+	    	"CCanvas.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(CCanvasError, 
+	    		"CCanvas.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 = CCanvas_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CCanvas_readDouble(PyList_GET_ITEM(pars,1));
+	    // Skip zero
+	    if (q_value==0) {
+	    	return Py_BuildValue("d",0.0);
+	    }
+		return Py_BuildValue("d",canvas_intensity(&(self->canvas_pars),q_value,phi_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CCanvas_readDouble(pars);		
+		
+		return Py_BuildValue("d",0.0);
+	}	
+}
+
+static PyObject * add(CCanvas *self, PyObject *args) {
+	int id;
+	
+	id = canvas_add(&(self->canvas_pars),REALSPACE_SPHERE);
+    return Py_BuildValue("i",id);
+}
+
+static PyObject * setParam(CCanvas *self, PyObject *args) {
+	int id, par_id;
+	double value;
+	double radius;
+
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"iid",&id, &par_id, &value) ) {
+	    PyErr_SetString(CCanvasError, 
+	    	"CCanvas.run expects a q value.");
+		return NULL;
+	}
+
+	canvas_setParam(&(self->canvas_pars),id,par_id,value);
+	
+	return Py_BuildValue("i",id);
+}
+
+
+static PyMethodDef CCanvas_methods[] = {
+    {"evaluate",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q"},
+    {"add",    (PyCFunction)add   , METH_VARARGS, "Add a shape object"},
+    {"setParam",    (PyCFunction)setParam   , METH_VARARGS, "Set a parameter"},
+   {NULL}
+};
+
+static PyTypeObject CCanvasType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CCanvas",             /*tp_name*/
+    sizeof(CCanvas),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CCanvas_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*/
+    "CCanvas objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CCanvas_methods,             /* tp_methods */
+    CCanvas_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CCanvas_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CCanvas_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 addCCanvas(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CCanvasType) < 0)
+        return;
+
+    Py_INCREF(&CCanvasType);
+    PyModule_AddObject(module, "CCanvas", (PyObject *)&CCanvasType);
+    
+    d = PyModule_GetDict(module);
+    CCanvasError = PyErr_NewException("CCanvas.error", NULL, NULL);
+    PyDict_SetItemString(d, "CCanvasError", CCanvasError);
+}
+
Index: sansmodels/src/sans/models/prototypes/src/CDispCylinderModel.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/CDispCylinderModel.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/CDispCylinderModel.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,256 @@
+/** CDispCylinderModel
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY disperse_cylinder.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 "disperse_cylinder.h"
+
+static double UNDEFINED = 1000010010.1234567;
+
+/// Error object for raised exceptions
+static PyObject * CDispCylinderModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Log for unit testing
+    PyObject * log;
+    /// Model parameters
+	DispCylinderParameters model_pars;
+} CDispCylinderModel;
+
+
+static void
+CDispCylinderModel_dealloc(CDispCylinderModel* self)
+{
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CDispCylinderModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CDispCylinderModel *self;
+    
+    self = (CDispCylinderModel *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CDispCylinderModel_init(CDispCylinderModel *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,"sigma_phi",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"sigma_radius",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"sigma_theta",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.000000));
+        PyDict_SetItemString(self->params,"cyl_theta",Py_BuildValue("d",1.570000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.000000));
+        PyDict_SetItemString(self->params,"n_pts",Py_BuildValue("d",25.000000));
+        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000003));
+        PyDict_SetItemString(self->params,"cyl_phi",Py_BuildValue("d",0.000000));
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CDispCylinderModel_members[] = {
+    {"params", T_OBJECT, offsetof(CDispCylinderModel, params), 0,
+     "Parameters"},
+    {"log", T_OBJECT, offsetof(CDispCylinderModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CDispCylinderModel_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(CDispCylinderModel *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.sigma_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma_phi") );
+    self->model_pars.sigma_radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma_radius") );
+    self->model_pars.sigma_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma_theta") );
+    self->model_pars.length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model_pars.cyl_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_theta") );
+    self->model_pars.background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model_pars.radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model_pars.n_pts = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_pts") );
+    self->model_pars.contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") );
+    self->model_pars.cyl_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_phi") );
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CDispCylinderModelError, 
+	    	"CDispCylinderModel.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(CDispCylinderModelError, 
+	    		"CDispCylinderModel.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 = CDispCylinderModel_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CDispCylinderModel_readDouble(PyList_GET_ITEM(pars,1));
+	    // Skip zero
+	    if (q_value==0) {
+	    	return Py_BuildValue("d",0.0);
+	    }
+		return Py_BuildValue("d",disperse_cylinder_analytical_2D(&(self->model_pars),q_value,phi_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CDispCylinderModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",disperse_cylinder_analytical_1D(&(self->model_pars),q_value));
+	}	
+}
+
+static PyObject * reset(CDispCylinderModel *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+
+static PyMethodDef CDispCylinderModel_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q"},
+    {"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 CDispCylinderModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CDispCylinderModel",             /*tp_name*/
+    sizeof(CDispCylinderModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CDispCylinderModel_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*/
+    "CDispCylinderModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CDispCylinderModel_methods,             /* tp_methods */
+    CDispCylinderModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CDispCylinderModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CDispCylinderModel_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 addCDispCylinderModel(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CDispCylinderModelType) < 0)
+        return;
+
+    Py_INCREF(&CDispCylinderModelType);
+    PyModule_AddObject(module, "CDispCylinderModel", (PyObject *)&CDispCylinderModelType);
+    
+    d = PyModule_GetDict(module);
+    CDispCylinderModelError = PyErr_NewException("CDispCylinderModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CDispCylinderModelError", CDispCylinderModelError);
+}
+
Index: sansmodels/src/sans/models/prototypes/src/CSimCylinder.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/CSimCylinder.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/CSimCylinder.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,246 @@
+/** CSimCylinder
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY simcylinder.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 "simcylinder.h"
+
+static double UNDEFINED = 1000010010.1234567;
+
+/// Error object for raised exceptions
+static PyObject * CSimCylinderError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Log for unit testing
+    PyObject * log;
+    /// Model parameters
+	SimCylinderParameters model_pars;
+} CSimCylinder;
+
+
+static void
+CSimCylinder_dealloc(CSimCylinder* self)
+{
+        modelcalculations_dealloc(&(self->model_pars.calcPars));
+    self->ob_type->tp_free((PyObject*)self);
+
+}
+
+static PyObject *
+CSimCylinder_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CSimCylinder *self;
+    
+    self = (CSimCylinder *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CSimCylinder_init(CSimCylinder *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"phi",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"qmax",Py_BuildValue("d",0.200000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",200.000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",80.000000));
+        PyDict_SetItemString(self->params,"theta",Py_BuildValue("d",1.570000));
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+                modelcalculations_init(&(self->model_pars.calcPars));
+
+    }
+    return 0;
+}
+
+static PyMemberDef CSimCylinder_members[] = {
+    {"params", T_OBJECT, offsetof(CSimCylinder, params), 0,
+     "Parameters"},
+    {"log", T_OBJECT, offsetof(CSimCylinder, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CSimCylinder_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(CSimCylinder *self, PyObject *args) {
+	double q_value, phi_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	// Reader parameter dictionary
+    self->model_pars.phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "phi") );
+    self->model_pars.scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model_pars.qmax = PyFloat_AsDouble( PyDict_GetItemString(self->params, "qmax") );
+    self->model_pars.length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model_pars.radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model_pars.theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "theta") );
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CSimCylinderError, 
+	    	"CSimCylinder.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(CSimCylinderError, 
+	    		"CSimCylinder.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 = CSimCylinder_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CSimCylinder_readDouble(PyList_GET_ITEM(pars,1));
+	    // Skip zero
+	    if (q_value==0) {
+	    	return Py_BuildValue("d",0.0);
+	    }
+		return Py_BuildValue("d",simcylinder_analytical_2D(&(self->model_pars),q_value,phi_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CSimCylinder_readDouble(pars);		
+		
+		return Py_BuildValue("d",simcylinder_analytical_1D(&(self->model_pars),q_value));
+	}	
+}
+
+static PyObject * reset(CSimCylinder *self, PyObject *args) {
+    modelcalculations_reset(&(self->model_pars.calcPars));
+
+    return Py_BuildValue("d",0.0);
+}
+
+
+static PyMethodDef CSimCylinder_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q"},
+    {"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 CSimCylinderType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CSimCylinder",             /*tp_name*/
+    sizeof(CSimCylinder),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CSimCylinder_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*/
+    "CSimCylinder objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CSimCylinder_methods,             /* tp_methods */
+    CSimCylinder_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CSimCylinder_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CSimCylinder_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 addCSimCylinder(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CSimCylinderType) < 0)
+        return;
+
+    Py_INCREF(&CSimCylinderType);
+    PyModule_AddObject(module, "CSimCylinder", (PyObject *)&CSimCylinderType);
+    
+    d = PyModule_GetDict(module);
+    CSimCylinderError = PyErr_NewException("CSimCylinder.error", NULL, NULL);
+    PyDict_SetItemString(d, "CSimCylinderError", CSimCylinderError);
+}
+
Index: sansmodels/src/sans/models/prototypes/src/CSimCylinderF.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/CSimCylinderF.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/CSimCylinderF.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,246 @@
+/** CSimCylinderF
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY ..\prototypes\src\simcylinder_fast.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 "simcylinder_fast.h"
+
+/// Error object for raised exceptions
+static PyObject * CSimCylinderFError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Log for unit testing
+    PyObject * log;
+    /// Model parameters
+	SimCylinderFParameters model_pars;
+} CSimCylinderF;
+
+
+static void
+CSimCylinderF_dealloc(CSimCylinderF* self)
+{
+    self->ob_type->tp_free((PyObject*)self);
+        modelcalculations_dealloc(&(self->model_pars.calcPars));
+
+}
+
+static PyObject *
+CSimCylinderF_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CSimCylinderF *self;
+    
+    self = (CSimCylinderF *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CSimCylinderF_init(CSimCylinderF *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	
+    	// Create parameters
+        self->params = PyDict_New();
+        
+        // Initialize parameter dictionary
+        PyDict_SetItemString(self->params,"theta",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.000000));
+        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.000000));
+        PyDict_SetItemString(self->params,"phi",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"npoints",Py_BuildValue("d",50000.000000));
+        PyDict_SetItemString(self->params,"seed",Py_BuildValue("d",1000.000000));
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+                modelcalculations_init(&(self->model_pars.calcPars));
+
+    }
+    return 0;
+}
+
+static PyMemberDef CSimCylinderF_members[] = {
+    {"params", T_OBJECT, offsetof(CSimCylinderF, params), 0,
+     "Parameters"},
+    {"log", T_OBJECT, offsetof(CSimCylinderF, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CSimCylinderF_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(CSimCylinderF *self, PyObject *args) {
+	double q_value, phi_value;
+	PyObject* pars;
+	int npars;
+	
+	// Get parameters
+	
+	// Reader parameter dictionary
+    self->model_pars.theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "theta") );
+    self->model_pars.length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model_pars.scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
+    self->model_pars.radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model_pars.phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "phi") );
+    self->model_pars.npoints = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npoints") );
+    self->model_pars.seed = PyFloat_AsDouble( PyDict_GetItemString(self->params, "seed") );
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CSimCylinderFError, 
+	    	"CSimCylinderF.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(CSimCylinderFError, 
+	    		"CSimCylinderF.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 = CSimCylinderF_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CSimCylinderF_readDouble(PyList_GET_ITEM(pars,1));
+	    // Skip zero
+	    if (q_value==0) {
+	    	return Py_BuildValue("d",0.0);
+	    }
+		return Py_BuildValue("d",simcylinder_fast_analytical_2D(&(self->model_pars),q_value,phi_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CSimCylinderF_readDouble(pars);		
+		
+		return Py_BuildValue("d",simcylinder_fast_analytical_1D(&(self->model_pars),q_value));
+	}	
+}
+
+static PyObject * reset(CSimCylinderF *self, PyObject *args) {
+    modelcalculations_reset(&(self->model_pars.calcPars));
+
+    return Py_BuildValue("d",0.0);
+}
+
+
+static PyMethodDef CSimCylinderF_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q"},
+    {"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 CSimCylinderFType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CSimCylinderF",             /*tp_name*/
+    sizeof(CSimCylinderF),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CSimCylinderF_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*/
+    "CSimCylinderF objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CSimCylinderF_methods,             /* tp_methods */
+    CSimCylinderF_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CSimCylinderF_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CSimCylinderF_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 addCSimCylinderF(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CSimCylinderFType) < 0)
+        return;
+
+    Py_INCREF(&CSimCylinderFType);
+    PyModule_AddObject(module, "CSimCylinderF", (PyObject *)&CSimCylinderFType);
+    
+    d = PyModule_GetDict(module);
+    CSimCylinderFError = PyErr_NewException("CSimCylinderF.error", NULL, NULL);
+    PyDict_SetItemString(d, "CSimCylinderFError", CSimCylinderFError);
+}
+
Index: sansmodels/src/sans/models/prototypes/src/CSimSphereF.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/CSimSphereF.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/CSimSphereF.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,240 @@
+/** CSimSphereF
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY ../prototypes/src/sphere_fast.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 "sphere_fast.h"
+
+/// Error object for raised exceptions
+static PyObject * CSimSphereFError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Log for unit testing
+    PyObject * log;
+    /// Model parameters
+	SimSphereFParameters model_pars;
+} CSimSphereF;
+
+
+static void
+CSimSphereF_dealloc(CSimSphereF* self)
+{
+    self->ob_type->tp_free((PyObject*)self);
+        modelcalculations_dealloc(&(self->model_pars.calcPars));
+
+}
+
+static PyObject *
+CSimSphereF_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CSimSphereF *self;
+    
+    self = (CSimSphereF *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CSimSphereF_init(CSimSphereF *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,"radius",Py_BuildValue("d",20.000000));
+        PyDict_SetItemString(self->params,"npoints",Py_BuildValue("d",50000.000000));
+        PyDict_SetItemString(self->params,"seed",Py_BuildValue("d",1000.000000));
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+                modelcalculations_init(&(self->model_pars.calcPars));
+
+    }
+    return 0;
+}
+
+static PyMemberDef CSimSphereF_members[] = {
+    {"params", T_OBJECT, offsetof(CSimSphereF, params), 0,
+     "Parameters"},
+    {"log", T_OBJECT, offsetof(CSimSphereF, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CSimSphereF_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(CSimSphereF *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.radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model_pars.npoints = PyFloat_AsDouble( PyDict_GetItemString(self->params, "npoints") );
+    self->model_pars.seed = PyFloat_AsDouble( PyDict_GetItemString(self->params, "seed") );
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CSimSphereFError, 
+	    	"CSimSphereF.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(CSimSphereFError, 
+	    		"CSimSphereF.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 = CSimSphereF_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CSimSphereF_readDouble(PyList_GET_ITEM(pars,1));
+	    // Skip zero
+	    if (q_value==0) {
+	    	return Py_BuildValue("d",0.0);
+	    }
+		return Py_BuildValue("d",sphere_fast_analytical_2D(&(self->model_pars),q_value,phi_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CSimSphereF_readDouble(pars);		
+		
+		return Py_BuildValue("d",sphere_fast_analytical_1D(&(self->model_pars),q_value));
+	}	
+}
+
+static PyObject * reset(CSimSphereF *self, PyObject *args) {
+    modelcalculations_reset(&(self->model_pars.calcPars));
+
+    return Py_BuildValue("d",0.0);
+}
+
+
+static PyMethodDef CSimSphereF_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q"},
+    {"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 CSimSphereFType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CSimSphereF",             /*tp_name*/
+    sizeof(CSimSphereF),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CSimSphereF_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*/
+    "CSimSphereF objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CSimSphereF_methods,             /* tp_methods */
+    CSimSphereF_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CSimSphereF_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CSimSphereF_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 addCSimSphereF(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CSimSphereFType) < 0)
+        return;
+
+    Py_INCREF(&CSimSphereFType);
+    PyModule_AddObject(module, "CSimSphereF", (PyObject *)&CSimSphereFType);
+    
+    d = PyModule_GetDict(module);
+    CSimSphereFError = PyErr_NewException("CSimSphereF.error", NULL, NULL);
+    PyDict_SetItemString(d, "CSimSphereFError", CSimSphereFError);
+}
+
Index: sansmodels/src/sans/models/prototypes/src/CSmearCylinderModel.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/CSmearCylinderModel.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/CSmearCylinderModel.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,254 @@
+/** CSmearCylinderModel
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY smeared_cylinder.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 "smeared_cylinder.h"
+
+static double UNDEFINED = 1000010010.1234567;
+
+/// Error object for raised exceptions
+static PyObject * CSmearCylinderModelError = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Log for unit testing
+    PyObject * log;
+    /// Model parameters
+	SmearCylinderParameters model_pars;
+} CSmearCylinderModel;
+
+
+static void
+CSmearCylinderModel_dealloc(CSmearCylinderModel* self)
+{
+    self->ob_type->tp_free((PyObject*)self);
+    
+
+}
+
+static PyObject *
+CSmearCylinderModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CSmearCylinderModel *self;
+    
+    self = (CSmearCylinderModel *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CSmearCylinderModel_init(CSmearCylinderModel *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,"sigma_phi",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"sigma_radius",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"sigma_theta",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"length",Py_BuildValue("d",400.000000));
+        PyDict_SetItemString(self->params,"cyl_theta",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.000000));
+        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000003));
+        PyDict_SetItemString(self->params,"cyl_phi",Py_BuildValue("d",0.000000));
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+        
+
+    }
+    return 0;
+}
+
+static PyMemberDef CSmearCylinderModel_members[] = {
+    {"params", T_OBJECT, offsetof(CSmearCylinderModel, params), 0,
+     "Parameters"},
+    {"log", T_OBJECT, offsetof(CSmearCylinderModel, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CSmearCylinderModel_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(CSmearCylinderModel *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.sigma_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma_phi") );
+    self->model_pars.sigma_radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma_radius") );
+    self->model_pars.sigma_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma_theta") );
+    self->model_pars.length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") );
+    self->model_pars.cyl_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_theta") );
+    self->model_pars.background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
+    self->model_pars.radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model_pars.contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") );
+    self->model_pars.cyl_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_phi") );
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CSmearCylinderModelError, 
+	    	"CSmearCylinderModel.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(CSmearCylinderModelError, 
+	    		"CSmearCylinderModel.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 = CSmearCylinderModel_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CSmearCylinderModel_readDouble(PyList_GET_ITEM(pars,1));
+	    // Skip zero
+	    if (q_value==0) {
+	    	return Py_BuildValue("d",0.0);
+	    }
+		return Py_BuildValue("d",smeared_cylinder_analytical_2D(&(self->model_pars),q_value,phi_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CSmearCylinderModel_readDouble(pars);		
+		
+		return Py_BuildValue("d",smeared_cylinder_analytical_1D(&(self->model_pars),q_value));
+	}	
+}
+
+static PyObject * reset(CSmearCylinderModel *self, PyObject *args) {
+    
+
+    return Py_BuildValue("d",0.0);
+}
+
+
+static PyMethodDef CSmearCylinderModel_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q"},
+    {"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 CSmearCylinderModelType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CSmearCylinderModel",             /*tp_name*/
+    sizeof(CSmearCylinderModel),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CSmearCylinderModel_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*/
+    "CSmearCylinderModel objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CSmearCylinderModel_methods,             /* tp_methods */
+    CSmearCylinderModel_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CSmearCylinderModel_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CSmearCylinderModel_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 addCSmearCylinderModel(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CSmearCylinderModelType) < 0)
+        return;
+
+    Py_INCREF(&CSmearCylinderModelType);
+    PyModule_AddObject(module, "CSmearCylinderModel", (PyObject *)&CSmearCylinderModelType);
+    
+    d = PyModule_GetDict(module);
+    CSmearCylinderModelError = PyErr_NewException("CSmearCylinderModel.error", NULL, NULL);
+    PyDict_SetItemString(d, "CSmearCylinderModelError", CSmearCylinderModelError);
+}
+
Index: sansmodels/src/sans/models/prototypes/src/CTestSphere2.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/CTestSphere2.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/CTestSphere2.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,240 @@
+/** CTestSphere2
+ *
+ * C extension 
+ *
+ * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
+ *          DO NOT MODIFY THIS FILE, MODIFY testsphere.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 "testsphere.h"
+
+static double UNDEFINED = 1000010010.1234567;
+
+/// Error object for raised exceptions
+static PyObject * CTestSphere2Error = NULL;
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD
+    /// Parameters
+    PyObject * params;
+    /// Log for unit testing
+    PyObject * log;
+    /// Model parameters
+	TestSphereParameters model_pars;
+} CTestSphere2;
+
+
+static void
+CTestSphere2_dealloc(CTestSphere2* self)
+{
+    self->ob_type->tp_free((PyObject*)self);
+        modelcalculations_dealloc(&(self->model_pars.calcPars));
+
+}
+
+static PyObject *
+CTestSphere2_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    CTestSphere2 *self;
+    
+    self = (CTestSphere2 *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+CTestSphere2_init(CTestSphere2 *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,"radius",Py_BuildValue("d",20.000000));
+        PyDict_SetItemString(self->params,"qmax",Py_BuildValue("d",0.100000));
+
+         
+        // Create empty log
+        self->log = PyDict_New();
+        
+                modelcalculations_init(&(self->model_pars.calcPars));
+
+    }
+    return 0;
+}
+
+static PyMemberDef CTestSphere2_members[] = {
+    {"params", T_OBJECT, offsetof(CTestSphere2, params), 0,
+     "Parameters"},
+    {"log", T_OBJECT, offsetof(CTestSphere2, log), 0,
+     "Log"},
+    {NULL}  /* Sentinel */
+};
+
+/** Read double from PyObject
+    @param p PyObject
+    @return double
+*/
+double CTestSphere2_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(CTestSphere2 *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.radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
+    self->model_pars.qmax = PyFloat_AsDouble( PyDict_GetItemString(self->params, "qmax") );
+
+	
+	// Get input and determine whether we have to supply a 1D or 2D return value.
+	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
+	    PyErr_SetString(CTestSphere2Error, 
+	    	"CTestSphere2.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(CTestSphere2Error, 
+	    		"CTestSphere2.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 = CTestSphere2_readDouble(PyList_GET_ITEM(pars,0));
+	    phi_value = CTestSphere2_readDouble(PyList_GET_ITEM(pars,1));
+	    // Skip zero
+	    if (q_value==0) {
+	    	return Py_BuildValue("d",0.0);
+	    }
+		return Py_BuildValue("d",testsphere_analytical_2D(&(self->model_pars),q_value,phi_value));
+
+	} else {
+
+		// We have a scalar q, we will evaluate I(q)
+		q_value = CTestSphere2_readDouble(pars);		
+		
+		return Py_BuildValue("d",testsphere_analytical_1D(&(self->model_pars),q_value));
+	}	
+}
+
+static PyObject * reset(CTestSphere2 *self, PyObject *args) {
+    modelcalculations_reset(&(self->model_pars.calcPars));
+
+    return Py_BuildValue("d",0.0);
+}
+
+
+static PyMethodDef CTestSphere2_methods[] = {
+    {"run",      (PyCFunction)run     , METH_VARARGS,
+      "Evaluate the model at a given Q"},
+    {"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 CTestSphere2Type = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "CTestSphere2",             /*tp_name*/
+    sizeof(CTestSphere2),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)CTestSphere2_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*/
+    "CTestSphere2 objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    CTestSphere2_methods,             /* tp_methods */
+    CTestSphere2_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)CTestSphere2_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    CTestSphere2_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 addCTestSphere2(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CTestSphere2Type) < 0)
+        return;
+
+    Py_INCREF(&CTestSphere2Type);
+    PyModule_AddObject(module, "CTestSphere2", (PyObject *)&CTestSphere2Type);
+    
+    d = PyModule_GetDict(module);
+    CTestSphere2Error = PyErr_NewException("CTestSphere2.error", NULL, NULL);
+    PyDict_SetItemString(d, "CTestSphere2Error", CTestSphere2Error);
+}
+
Index: sansmodels/src/sans/models/prototypes/src/c_models.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/c_models.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/c_models.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,39 @@
+/** c_models
+ *
+ * Module containing all SANS model extensions 
+ *
+ * @author   M.Doucet / UTK
+ */
+#include <Python.h>
+
+/// CVS information
+static char cvsid[] = "$Id: BaseCComponent.c,v 1.4 2007/03/23 20:50:58 doucet Exp $";
+
+
+/**
+ * Define empty module
+ */
+static PyMethodDef module_methods[] = {
+    {NULL} 
+};
+
+
+#ifndef PyMODINIT_FUNC	/* declarations for DLL import/export */
+#define PyMODINIT_FUNC void
+#endif
+PyMODINIT_FUNC
+initc_models(void) 
+{
+    PyObject* m;
+
+    m = Py_InitModule3("c_models", module_methods,
+                       "C extension module for SANS scattering models.");
+                       
+	addCTestSphere2(m);
+	addCSimCylinder(m);
+	addCSmearCylinderModel(m);
+	addCDispCylinderModel(m);
+	addCSimCylinderF(m);
+	addCSimSphereF(m);
+	addCCanvas(m);
+}
Index: sansmodels/src/sans/models/prototypes/src/canvas.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/canvas.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/canvas.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,291 @@
+#include "modelCalculations.h"
+#include "canvas.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <memory.h>
+#include <time.h>
+
+/**
+ * Initialization function for simulation structure
+ */
+int canvas_init(CanvasParams *pars) {
+	int i;
+	
+	pars->n_shapes = 0;
+	pars->max_shapes = 3; 
+	
+	// Allocate memory
+	pars->shapes = (SpaceObject*) malloc(pars->max_shapes*sizeof(SpaceObject));
+	
+	if(pars->shapes==NULL) {
+		printf("canvas_init: Problem allocating memory for space objects\n");
+		return -1;
+	}
+	
+	return 0;
+}
+
+int canvas_dealloc(CanvasParams *self) {
+	int i;
+	printf("dealloc\n");
+	/*
+	//for(i=0; i<self->n_shapes; i++) {
+	for(i=0; i<self->max_shapes; i++) {
+		free(self->shapes[i].points);
+		free(self->shapes[i].params);
+		//free(self->shapes[i].generate);
+	}
+	*/
+	free(self->shapes);
+	printf("deallocated\n");
+	return 0;
+}
+
+int canvas_add(CanvasParams *self, int objectCode) {
+	int id;
+	id = self->n_shapes;
+	
+	if(objectCode==REALSPACE_SPHERE) {
+		initSphereObject( &(self->shapes[id]) );
+	}
+	
+	self->shapes[id].layer = id;
+	self->n_shapes++;
+	
+	return id;
+}
+
+int canvas_setParam(CanvasParams *self, int shapeID, int paramID, double value) {
+	self->shapes[shapeID].params[paramID]=value;
+	self->shapes[shapeID].points_available = 0;
+	
+	// update volume
+	self->shapes[shapeID].volume = (*self->shapes[shapeID].getVolume)(&(self->shapes[shapeID].params));
+	
+	return 0;
+}
+
+double canvas_intensity(CanvasParams *self, double q, double phi) {
+	double phase, cos_term, sin_term;
+	double qx, qy, vol;
+	int i, j, npts;
+	double pars[4];
+	
+	double cos_shape, sin_shape;
+	
+	double vol_core, vol_shell, vol_sum;
+	
+	pars[0] = 1.0;
+	pars[1] = 40.0;
+	
+	cos_term = 0.0;
+	sin_term = 0.0;
+	
+	vol_sum = 0.0;
+	
+	vol = self->shapes[0].volume;
+	vol_core = self->shapes[1].volume;
+	vol_shell = vol-vol_core;
+
+	npts = 0;
+	
+	qx = q*cos(phi);
+	qy = q*sin(phi);
+	
+	for(i=0; i<self->n_shapes; i++) {
+		if (self->shapes[i].points_available==0) {
+			(*self->shapes[i].generate)( self->shapes[i].points, self->shapes[i].params, self->shapes[i].npts);
+			self->shapes[i].points_available=1;
+		}
+		npts = 0;
+		cos_shape = 0;
+		sin_shape = 0;
+		for(j=0; j<self->shapes[i].npts; j++) {
+			
+			if(canvas_PointAllowed(self, i, j)==1) {
+				npts++;
+				
+				phase = qx*self->shapes[i].points[j].x + qy*self->shapes[i].points[j].y;
+				//phase = 0.0;
+				
+				cos_shape += cos(phase);
+				sin_shape += sin(phase);
+			}   			
+		}
+
+		//printf("ID %i: cos=%g sin=%g fac=%g\n",i,cos_shape, sin_shape, self->shapes[i].volume /self->shapes[i].npts  * self->shapes[i].params[2]);
+
+
+		
+		
+		// TODO: must find an efficient way to get the exact volume of each SLD regions.
+		// For instance, if a volume is completely included in another (and it has a higher layer
+		// number, we can simply subtract the two "theoretical" volumes instead of using the 
+		// ratio of the points.
+
+		// The following are the three lines we would like to write:
+		
+		//cos_term += cos_shape * self->shapes[i].volume /npts  * self->shapes[i].params[2];
+		//sin_term += sin_shape * self->shapes[i].volume /npts  * self->shapes[i].params[2];
+		//vol_sum += self->shapes[i].volume;
+		
+		// In that case, the volume is approx
+		//     vol = (full volume) - sum(volume parts overlapping with other objects)
+		// That can easily be done for object completely contained by others.
+
+		cos_term += cos_shape * self->shapes[i].volume /self->shapes[i].npts  * self->shapes[i].params[2];
+		sin_term += sin_shape * self->shapes[i].volume /self->shapes[i].npts  * self->shapes[i].params[2];
+		vol_sum += self->shapes[i].volume*npts/self->shapes[i].npts;
+		
+		
+		//vol_sum += self->shapes[i].volume;
+
+		// The following is more precise (but implies better investigation of the topology)
+		/*
+		if(i==0){
+			cos_term += cos_shape * vol_shell/npts  * self->shapes[i].params[2];
+			sin_term += sin_shape * vol_shell/npts  * self->shapes[i].params[2];
+			vol_sum += vol_shell;
+		} else {
+			cos_term += cos_shape * vol_core/npts  * self->shapes[i].params[2];
+			sin_term += sin_shape * vol_core/npts  * self->shapes[i].params[2];
+			vol_sum += vol_core;
+		}
+		*/
+	}
+	
+	return 1.0e8*(cos_term*cos_term + sin_term*sin_term)/vol_sum;	
+} 
+
+int canvas_PointAllowed(CanvasParams *self, int object_id, int i_pt) {
+	// Check whether space point i_pt is already covered by another object
+	int i;
+	
+	for(i=0; i<self->n_shapes; i++) {
+		// Check the object object_id is underneath
+		if(self->shapes[object_id].layer < self->shapes[i].layer) {
+			// Is the point overlapping, if so skip it.
+			if(self->shapes[i].isContained(&self->shapes[object_id].points[i_pt], 
+				self->shapes[i].params) == 1) {
+				return 0;
+			}
+		} 
+	}
+	return 1;
+} 
+
+
+
+// SPHERE object ----------------------------------------------------------------------
+
+
+int initSphereObject(SpaceObject *pars) {
+	
+	pars->x = 0.0;
+	pars->y = 0.0;
+	pars->z = 0.0;
+	
+	pars->o1 = 0.0;
+	pars->o2 = 0.0;
+	pars->o3 = 0.0;
+	
+	/*
+	free(pars->params);
+	
+	printf("allocating parameters\n");
+	pars->params = (double*) malloc(5*sizeof(double));
+	if(pars->params==NULL) {
+		printf("initSphereObject: Problem allocating memory\n");
+		return -1;	
+	}
+	*/
+
+	// Scale
+	pars->params[0] = 1.0;
+	// Radius
+	pars->params[1] = 40.0;
+	// Contrast
+	pars->params[2] = 1.0;
+	 
+	pars->points_available = 0;
+	
+	pars->npts = 150000;
+	free(pars->points);
+	
+	// Generation function
+	pars->generate = &generateSphereObject;
+	// isContained function
+	pars->isContained = &sphere_IsContained;
+	pars->getVolume = &sphere_getVolume;
+	
+	pars->volume = acos(-1.0)*4.0/3.0*pars->params[1]*pars->params[1]*pars->params[1];
+	
+	pars->layer = 0;
+	
+	return 0;
+	
+}
+
+int sphere_IsContained(SpacePoint *point, double *params) {
+	if( sqrt(point->x*point->x + point->y*point->y + point->z*point->z) < params[1]) {
+	  return 1;
+	}
+	return 0;
+}
+
+double sphere_getVolume(double *params) {
+	printf("radius=%g\n",params[1]);
+	return acos(-1.0)*4.0/3.0*params[1]*params[1]*params[1];
+}
+
+
+int generateSphereObject(SpacePoint *points, double *params, int npts) {
+	int i, testcounter;
+	double x, y, z;
+	SpacePoint tmp;
+	
+	printf("radius=%g, npts=%i\n", params[1], npts);
+	
+	testcounter = 0;
+	
+	/*
+	free(points);
+	printf("allocating points\n");
+	points = (SpacePoint*) malloc(npts*sizeof(SpacePoint));
+	if(points==NULL) {
+		printf("generateSphereObject: Problem allocating memory\n");
+		return -1;	
+	}
+	*/
+	memset(points,0,npts*sizeof(SpacePoint));
+
+	for(i=0;i<npts;i++) {
+		// Generate in a box centered around zero
+		x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * params[1];
+		y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * params[1];
+		z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * params[1];
+		
+		// reject those that are not within the volume
+//		if( sqrt(x*x+y*y+z*z) < params[1]) {
+		tmp.x = x;
+		tmp.y = y;
+		tmp.z = z;
+		if( sphere_IsContained(&tmp, params) == 1 ) {
+			points[i].x =  x;
+			points[i].y =  y;
+			points[i].z =  z;
+			testcounter++;
+		} else {
+			i--;
+		}
+	}
+	
+	// Consistency check
+	if(testcounter != npts) {
+		return -1;
+	}
+	
+	return testcounter;
+	
+}
Index: sansmodels/src/sans/models/prototypes/src/canvas.h
===================================================================
--- sansmodels/src/sans/models/prototypes/src/canvas.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/canvas.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,65 @@
+#if !defined(simcanvas_h)
+#define simcanvas_h
+
+#include "modelCalculations.h"
+
+#define REALSPACE_SPHERE 1
+
+typedef struct {
+	// Position
+	double x;
+	double y;
+	double z;
+	// Orientation
+	double o1;
+	double o2;
+	double o3;
+	// Vector of parameters
+	double params[15];
+	
+	// Pointer to function to generate space points
+	int (*generate)(SpacePoint *points, double *params, int npts);
+	// Pointer to function to check whether a point is contained by an object
+	int (*isContained)(SpacePoint *point, double *params);
+	// Pointer to function to calculate the object's volume
+	double (*getVolume)(double *params);
+	
+	
+	// Points available
+	int points_available;
+	// Space points
+	SpacePoint points[151000];
+	//SpacePoint *points;
+	// Number of space points
+	int npts;
+	double volume;
+	int layer;
+	
+} SpaceObject;
+
+
+typedef struct {
+	// Maximum number of shapes
+	int max_shapes;
+	// List of real-space objects
+	SpaceObject *shapes;
+	int n_shapes;
+} CanvasParams;
+
+
+int canvas_dealloc(CanvasParams *self);
+int canvas_add(CanvasParams *pars, SpaceObject *shape);
+double canvas_intensity(CanvasParams *self, double q, double phi);
+int canvas_init(CanvasParams *pars);
+int canvas_add(CanvasParams *self, int objectCode);
+int canvas_setParam(CanvasParams *self, int shapeID, int paramID, double value);
+int canvas_PointAllowed(CanvasParams *self, int object_id, int i_pt);
+
+
+// SPHERE object ----------------------------------------------------------------------
+int initSphereObject(SpaceObject *pars);
+int generateSphereObject(SpacePoint *points, double *params, int npts);
+int sphere_IsContained(SpacePoint *point, double *params);
+double sphere_getVolume(double *params);
+
+#endif
Index: sansmodels/src/sans/models/prototypes/src/disperse_cylinder.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/disperse_cylinder.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/disperse_cylinder.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,86 @@
+/**
+ * Scattering model for a cylinder
+ * @author: Mathieu Doucet / UTK
+ */
+
+#include "cylinder.h"
+#include "disperse_cylinder.h"
+#include <math.h>
+#include "libCylinder.h"
+
+double c_disperser( double (*eval)(), double dp[], int n_pars, 
+					int *idList, double *sigmaList, int n_pts, double q, double phi );
+
+/**
+ * Function to evaluate 1D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+double disperse_cylinder_analytical_1D(DispCylinderParameters *pars, double q) {
+	double dp[5];
+	
+	// Fill paramater array
+	dp[0] = pars->scale;
+	dp[1] = pars->radius;
+	dp[2] = pars->length;
+	dp[3] = pars->contrast;
+	dp[4] = pars->background;
+	
+	return CylinderForm(dp, q);
+}
+	
+double disp_cyl_run(double pars[], double q, double phi) {
+	CylinderParameters danse_pars;
+	danse_pars.scale        = pars[0];
+	danse_pars.radius       = pars[1];
+	danse_pars.length       = pars[2];
+	danse_pars.contrast     = pars[3];
+	danse_pars.background   = pars[4];
+	danse_pars.cyl_theta    = pars[5];
+	danse_pars.cyl_phi      = pars[6];
+
+	return cylinder_analytical_2D(&danse_pars, q, phi);		
+		
+}
+	
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+double disperse_cylinder_analytical_2D(DispCylinderParameters *pars, double q, double phi) {
+	CylinderParameters cyl_pars;
+	double dp[6];
+	int paramList[3];
+	double sigmaList[3];
+	int npts;
+	
+	
+	// Fill paramater struct
+	dp[0] = pars->scale;
+	dp[1] = pars->radius;
+	dp[2] = pars->length;
+	dp[3] = pars->contrast;
+	dp[4] = pars->background;
+	dp[6] = pars->cyl_phi;
+	dp[5] = pars->cyl_theta;
+	
+	paramList[0] = 6;
+	paramList[1] = 5;
+	paramList[2] = 1;
+	
+	
+	sigmaList[0] = pars->sigma_phi;
+	sigmaList[1] = pars->sigma_theta;
+	sigmaList[2] = pars->sigma_radius;
+	
+	npts = (int)(floor(pars->n_pts));
+	
+	return c_disperser( &disp_cyl_run, &dp, 3, 
+					paramList, sigmaList, npts, q, phi ); 
+
+	
+}
+    
Index: sansmodels/src/sans/models/prototypes/src/disperse_cylinder.h
===================================================================
--- sansmodels/src/sans/models/prototypes/src/disperse_cylinder.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/disperse_cylinder.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,67 @@
+#if !defined(disperse_cylinder_h)
+#define disperse_cylinder_h
+
+/** Structure definition for dispersed cylinder parameters 
+ * [PYTHONCLASS] = DispCylinderModel
+ * [BASESTRUCT]  = CylinderParameters 
+ * [BASE_H]      = cylinder
+ * */
+typedef struct {
+    /// Scale factor 
+    //  [DEFAULT]=scale=1.0
+    double scale;
+    /// Radius of the cylinder [A]
+    //  [DEFAULT]=radius=20.0 A
+    double radius;
+    /// Length of the cylinder [A]
+    //  [DEFAULT]=length=400.0 A
+    double length;
+    /// Contrast [A-2]
+    //  [DEFAULT]=contrast=3.0e-6 A-2
+    double contrast;
+	/// Incoherent Background (cm-1) 0.000
+	//  [DEFAULT]=background=0 cm-1
+	double background;    
+    /// Orientation of the cylinder axis w/respect incoming beam [rad]
+    //  [DEFAULT]=cyl_theta=1.57 rad
+    double cyl_theta;
+    /// Orientation of the cylinder in the plane of the detector [rad]
+    //  [DEFAULT]=cyl_phi=0.0 rad
+    
+    ///
+    //  [DEFAULT]=cyl_phi=0.0 rad
+    double cyl_phi;    
+    
+    ///
+    //  [DEFAULT]=sigma_theta=0.0 rad
+    double sigma_theta;
+    
+    ///
+    //  [DEFAULT]=sigma_phi=0.0 rad
+    double sigma_phi;
+    
+    ///
+    //  [DEFAULT]=sigma_radius=0.0 rad
+    double sigma_radius;
+    
+    ///
+    //  [DEFAULT]=sigma_radius=0.0 rad
+    double sigma_length;
+    
+    ///
+    //  [DEFAULT]=n_pts=25
+    double n_pts;
+    
+    
+} DispCylinderParameters;
+
+
+
+/// 1D scattering function
+double disperse_cylinder_analytical_1D(DispCylinderParameters *pars, double q);
+
+/// 2D scattering function
+double disperse_cylinder_analytical_2D(DispCylinderParameters *pars, double q, double phi);
+double disperse_cylinder_dist( double x, double mean, double sigma );
+
+#endif
Index: sansmodels/src/sans/models/prototypes/src/modelCalculations.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/modelCalculations.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/modelCalculations.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,445 @@
+#include "modelCalculations.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <memory.h>
+#include <time.h>
+
+/**
+ * Initialization function for simulation structure
+ */
+void modelcalculations_init(CalcParameters *pars) {
+        pars->isRhoAvailable         = 0;
+        pars->isPointMemAllocated    = 0;
+        pars->isRhoAvailable_2D      = 0;
+        pars->isPointMemAllocated_2D = 0;
+       
+       	pars->volume_points = 0;
+       	pars->r_points = 0;
+       
+		pars->timePr_1D = 0;
+    	pars->timePr_2D = 0;
+    	pars->timeIq_1D = 0;
+    	pars->timeIq_2D = 0;
+    
+    	pars->errorOccured = 0;
+    	
+}
+
+/**
+ * Reset function for simulation structure
+ */
+void modelcalculations_reset(CalcParameters *pars) {
+       modelcalculations_dealloc(pars);
+       modelcalculations_init(pars);
+}
+
+/**
+ * Deallocate memory of simullation structure
+ */
+void modelcalculations_dealloc(CalcParameters *pars) {
+	free(pars->rho);
+	free(pars->points);
+	free(pars->rho_2D);
+	free(pars->points_2D);
+}
+
+/**
+ * Calculate pair correlation for 1D simulation
+ */
+int modelcalculations_calculatePairCorrelation_1D(SpacePoint * points, int volume_points, double * rho, int r_points, double bin_width) {
+	int i,j;
+	double dx,dy,dz,dist;
+	int i_bin;
+	//double bin_width;
+	double delta_t;
+	double average;
+	double closest;
+	time_t start_time;
+	clock_t start;
+	clock_t finish;
+	//struct tm *timeStruct;
+	
+	
+	// Allocate memory
+	/*
+	rho = (double*) malloc(r_points*sizeof(double));
+	if(rho==NULL){
+		printf("Problem allocating memory for 1D correlation points\n");
+		return -1;
+	}
+	*/
+	
+	// Clear vector
+	memset(rho,0,r_points*sizeof(double));
+
+	// R bin width
+    //bin_width = 2.0*size/r_points;
+    
+		time(&start_time);
+		start = clock();
+		
+		
+	average = 0;
+	for(i=0;i<volume_points-1;i++) {
+		closest = -1;
+		for(j=i+1;j<volume_points;j++) {
+            dx = (points[i].x-points[j].x);
+            dy = (points[i].y-points[j].y);
+            dz = (points[i].z-points[j].z);
+			dist = sqrt(dx*dx + dy*dy + dz*dz);
+			
+			if(closest<0 || dist<closest) {
+				closest = dist;
+			}
+			
+			//i_bin = (int)dist/bin_width;
+			i_bin = (int)floor(dist/bin_width);
+            
+            //rho[i_bin] = rho[i_bin] + 1.0/9000000.0;
+            if(i_bin >= r_points) {
+            	printf("problem! %i > %i\n", i_bin, r_points);
+            } else {
+            	rho[i_bin] = rho[i_bin] + 1.0;
+            }
+		}
+		average += closest;
+	}
+	average = average/(double)volume_points;
+	printf("average distance %f\n",average);
+	
+		finish = clock();
+        //time(&end_time);
+        //delta_t = difftime(end_time,start_time);
+        delta_t = ((double)(finish-start))/CLOCKS_PER_SEC;
+        printf("------------->PR calc time = %f\n", delta_t);
+    return 1;
+}
+
+/**
+ * Calculate I(q) for 1D simulation 
+ */
+double modelcalculations_calculateIq_1D(double * rho, int r_points, double r_step, double q) {
+	int i;
+	double value;
+	//double r_step;
+	//double vol;
+	double sum;
+	double qr;
+	clock_t start;
+	clock_t finish;
+	double delta_t;
+	
+	start = clock();
+	
+    //vol = 4.0*acos(-1.0)/3.0*radius*radius*radius;
+	//r_step = 2.0*radius/((double)(r_points));
+	
+	value = 0.0;
+	sum = 0.0;
+	for(i=1; i<r_points; i++) {
+		qr = q*r_step*(double)i;
+		value = value + rho[i] * sin(qr) / qr;
+		sum = sum + rho[i];
+	}
+	
+	
+	value = value/sum;
+	
+	finish = clock();
+    delta_t = ((double)(finish-start))/CLOCKS_PER_SEC;
+    //printf("------------->IQ calc time = %f\n", delta_t);
+	
+	
+	return value;
+}
+
+/**
+ * Calculate pair correlation function for 2D simulation using
+ * a 3D array to store P(r)
+ */
+int modelcalculations_calculatePairCorrelation_2D_3Darray(SpacePoint * points, int volume_points, float * rho, int r_points, double bin_width) {
+	int i,j;
+	int ix_bin, iy_bin, iz_bin;
+	
+	clock_t start;
+	clock_t finish;
+	
+	// Clear vector
+	memset(rho,0,r_points*r_points*r_points*sizeof(float));
+
+    start = clock();
+	for(i=0;i<volume_points-1;i++) {
+		for(j=i+1;j<volume_points;j++) {
+			
+            // Add entry to the matrix
+            ix_bin = (int)floor(fabs(points[i].x-points[j].x)/bin_width);
+            iy_bin = (int)floor(fabs(points[i].y-points[j].y)/bin_width);
+            iz_bin = (int)floor(fabs(points[i].z-points[j].z)/bin_width);
+            
+            if(ix_bin < r_points && iy_bin < r_points && iz_bin < r_points) {
+            	rho[(ix_bin*r_points+iy_bin)*r_points+iz_bin] += 1.0;
+            } else {
+            	printf("Bad point! %i %i %i\n", ix_bin, iy_bin, iz_bin);
+            }
+			
+		}
+	}
+	
+	finish = clock();
+    printf("-------------> Pair Correlation time = %f\n", ((double)(finish-start))/CLOCKS_PER_SEC);
+    return 0;
+}
+
+/**
+ * Calculate pair correlation function for 2D simulation by storing P(r) in a 2D array.
+ * Allows for rotation of object in space by specifying theta, phi, omage of the beam
+ */
+int modelcalculations_calculatePairCorrelation_2D_vector(SpacePoint * points, int volume_points, float * rho, 
+			int r_points, double bin_width, double theta_beam, double phi_beam, double omega_beam) {
+	int i,j;
+	int ix_bin, iy_bin;
+	SpacePoint p1, p2;
+	
+	clock_t start;
+	clock_t finish;
+	
+	// Clear vector
+	memset(rho,0,r_points*r_points*sizeof(float));
+	//printf("P(r) with theta=%g phi=%g\n", theta_beam, phi_beam);
+    start = clock();
+	for(i=0;i<volume_points-1;i++) {
+		// Rotate point
+		p1 = modelcalculations_rotate(points[i], theta_beam, phi_beam, omega_beam);
+		//printf("p = %g %g %g -> %g %g %g\n", points[i].x,points[i].y,points[i].z, p1.x,p1.y,p1.z);
+		for(j=i+1;j<volume_points;j++) {
+			// Rotate point
+			p2 = modelcalculations_rotate(points[j], theta_beam, phi_beam, omega_beam);
+			
+			
+			// Calculate distance in plane perpendicular to beam (z)
+            ix_bin = (int)floor(fabs(p1.x-p2.x)/bin_width);
+            iy_bin = (int)floor(fabs(p1.y-p2.y)/bin_width);
+            //iz_bin = (int)floor((points[i].z-points[j].z)/bin_width+r_points/2);
+            
+            rho[ix_bin*r_points+iy_bin] += 1.0;
+			
+		}
+	}
+	
+	finish = clock();
+    printf("-------------> 2D (v) Pair Correlation time = %f\n", ((double)(finish-start))/CLOCKS_PER_SEC);
+    return 1;
+}
+
+/**
+ * Calculate pair correlation function for 2D simulation by storing P(r) in a 2D array
+ */
+int modelcalculations_calculatePairCorrelation_2D(SpacePoint * points, int volume_points, float * rho, int r_points, double bin_width) {
+	int i,j;
+	int ix_bin, iy_bin;
+	
+	clock_t start;
+	clock_t finish;
+	
+	// Clear vector
+	memset(rho,0,r_points*r_points*sizeof(float));
+
+	//return 1;
+	
+    start = clock();
+	for(i=0;i<volume_points-1;i++) {
+		for(j=i+1;j<volume_points;j++) {
+			
+			// Calculate distance in plane perpendicular to beam (z)
+            ix_bin = (int)floor(fabs(points[i].x-points[j].x)/bin_width);
+            iy_bin = (int)floor(fabs(points[i].y-points[j].y)/bin_width);
+            //iz_bin = (int)floor(fabs(points[i].z-points[j].z)/bin_width);
+            
+            rho[ix_bin*r_points+iy_bin] += 1.0;
+            //rho[ix_bin*r_points+iy_bin] += fabs(points[i].z-points[j].z);
+		}
+	}
+	
+	finish = clock();
+    printf("-------------> Pair Correlation time = %f\n", ((double)(finish-start))/CLOCKS_PER_SEC);
+    
+    /*
+     * for(i=0;i<r_points;i++) {
+    	for(j=0;j<r_points;j++) {
+    		printf("Pr(%i, %i) = %g\n", i, j, rho[i*r_points+j]);
+    	}	
+    	
+    }
+    */
+	return 0;
+}
+
+/**
+ * Calculate I(q) for 2D simulation from 3D array
+ */
+double modelcalculations_calculateIq_2D_3Darray(float * rho, int r_points, double r_step, double q, double phi) {
+	//TODO: make rho an array of ints.
+	
+	int ix,iy,iz;
+	
+	// This should be a parameter
+	double lambda = 1.6;
+	double theta;
+	double value;
+	double sum;
+	// This should also be a parameter, about value of radius
+	double r_max = 1;
+	//double r_step =1;
+	
+	double c1;
+	double c2;
+	double c3;
+	int r2;
+	double f3;
+	double iz_c;
+	clock_t start; 
+	clock_t finish;
+	
+	theta = 2*asin(q*lambda/(4*acos(-1.0)));
+	
+	value = 0.0;
+	sum = 0.0;
+	//c1 = lambda*(cos(theta)-1)*r_step;
+	//c2 = lambda*sin(theta)*cos(phi)*r_step;
+	//c3 = lambda*sin(theta)*sin(phi)*r_step;
+	
+	c1 = 0;
+	c2 = q*cos(phi)*r_step;
+	c3 = q*sin(phi)*r_step;
+	
+	start = clock();
+	
+	// TODO: sum(rho) should be equal to the number of points^2 !
+	
+	for(ix=0; ix<r_points; ix++) {
+		for(iy=0; iy<r_points; iy++) {
+			r2 = ix*r_points*r_points+iy*r_points;
+			f3 = c2*(ix-r_points/2) + c3*(iy-r_points/2);
+			for(iz=0; iz<r_points; iz++) {
+				iz_c = (double)iz-r_points/2; 
+				value = value + rho[r2+iz] * cos( c1*iz_c + f3 );
+				sum = sum + rho[r2+iz];
+			}
+		}
+	}
+	finish = clock();
+    //printf("-------------> I(Q) time = %f, (%f %f %f)\n", ((double)(finish-start))/CLOCKS_PER_SEC, q, phi,value/sum);
+	
+	
+	value = value /sum;
+	return value;
+}
+
+/**
+ * Pair correlation function for a sphere
+ * 	@param r: distance value
+ */
+double pair_corr_sphere(double r) {
+	
+	return r*r*(1.0 - 0.75*r + r*r*r/16.0);	
+}
+
+/** 
+ * Calculate I(q) for 2D simulation from 2D array P(r)
+ */
+double modelcalculations_calculateIq_2D(float * rho, int r_points, double r_step, double q, double phi) {
+	//TODO: make rho an array of ints.
+	
+	int ix,iy;
+	
+	// This should be a parameter
+	double lambda = 1.0;
+	double value;
+	double sum;
+	// This should also be a parameter, about value of radius
+	double r_max = 1;
+	//double r_step =1;
+	
+	double c2;
+	double c3;
+	int ibin;
+	clock_t start; 
+	clock_t finish;
+	
+	//theta = 2*asin(q*lambda/(4*acos(-1.0)));
+	value = 0.0;
+	sum = 0.0;
+	c2 = q*cos(phi)*r_step;
+	c3 = q*sin(phi)*r_step;
+	//printf("phi = %g, c2=%g, c3=%g, q=%g, r_step=%g\n",phi, c2, c3,q,r_step);
+	
+	start = clock();
+	
+	for(ix=-r_points+1; ix<r_points; ix++) {
+		for(iy=-r_points+1; iy<r_points; iy++) {
+			
+			//value += rho[ix*r_points+iy] * cos( c2*((double)ix+0.5) + c3*((double)iy+0.5) );
+			//sum += rho[ix*r_points+iy];
+			
+			ibin = (int)(floor(sqrt(1.0*ix*ix)))*r_points+(int)(floor(sqrt(1.0*iy*iy)));
+			
+			if (ibin<r_points*r_points) {
+			
+				//value += rho[ibin] * cos( c2*((double)ix+0.5) + c3*((double)iy+0.5) );
+			
+				value += rho[ibin] * cos( c2*((double)ix) + c3*((double)iy) );
+				
+				sum += rho[ibin];
+			} else {
+				printf("Error computing IQ %i >= %i (%i %i)\n", ibin, r_points*r_points,ix, iy);
+			};
+			
+			
+			//dx = r_step*((double)ix+0.5);
+			//dy = r_step*((double)iy+0.5);
+			
+			// Only works for sphere of radius = 20
+			//f3 = pair_corr_sphere(sqrt(dx*dx+dy*dy)/20.0);
+			
+			//value += f3 * cos( c2*(dx) + c3*(dy) );
+			//sum += f3;		
+		}
+	}
+	finish = clock();
+	
+	value = value /sum;
+	return value;
+}
+
+
+/**
+ *  Rotation of a space point
+ */
+SpacePoint modelcalculations_rotate(SpacePoint p, double theta, double phi, double omega) {
+	SpacePoint new_point;
+	double x_1, x_2;
+	double y_1, y_2;
+	double z_1, z_2;
+	// P(r) assumes beam along z-axis. Rotate point accordingly
+	
+	
+	// Omega, around z-axis (doesn't change anything for cylindrical symmetry
+	x_1 = p.x*cos(omega) - p.y*sin(omega);
+	y_1 = p.x*sin(omega) + p.y*cos(omega);
+	z_1 = p.z;	
+	
+	// Theta, around y-axis
+	x_2 = x_1*cos(theta) + z_1*sin(theta);
+	y_2 = y_1;
+	z_2 = -x_1*sin(theta) + z_1*cos(theta);
+	
+	// Phi, around z-axis
+	new_point.x = x_2*cos(phi) - y_2*sin(phi);
+	new_point.y = x_2*sin(phi) + y_2*cos(phi);
+	new_point.z = z_2;
+	
+	return new_point;	
+}
+ 
+ 
Index: sansmodels/src/sans/models/prototypes/src/modelCalculations.h
===================================================================
--- sansmodels/src/sans/models/prototypes/src/modelCalculations.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/modelCalculations.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,64 @@
+#if !defined(modelcalcs_h)
+#define modelcalcs_h
+
+typedef struct {
+	double x;
+	double y;
+	double z;
+} SpacePoint;
+
+
+
+typedef struct {
+    /// Volume points
+    SpacePoint * points;
+    /// Pair correlation function
+    double * rho;
+    /// Is rho already generated?
+    int isRhoAvailable;
+    /// Memory allocation flag 
+    int isPointMemAllocated;
+    
+    /// Volume points
+    SpacePoint * points_2D;
+    /// Pair correlation function
+    float * rho_2D;
+    /// Is rho already generated?
+    int isRhoAvailable_2D;
+    /// Memory allocation flag
+    int isPointMemAllocated_2D;
+    
+    //TODO: add volume_points, r_points, step here...
+    
+    int volume_points;
+    int r_points;
+    
+    float timePr_1D;
+    float timePr_2D;
+    float timeIq_1D;
+    float timeIq_2D;
+    
+    int errorOccured;
+    
+    
+} CalcParameters;
+
+void modelcalculations_dealloc(CalcParameters *pars);
+void modelcalculations_init(CalcParameters *pars);
+void modelcalculations_reset(CalcParameters *pars);
+
+
+double modelcalculations_calculateIq_1D(double * rho, int r_points, double step, double q);
+int modelcalculations_calculatePairCorrelation_1D(SpacePoint * points, int volume_points, double * rho, int r_points, double size);
+
+double modelcalculations_calculateIq_2D(float * rho, int r_points, double step, double q, double phi);
+int modelcalculations_calculatePairCorrelation_2D(SpacePoint * points, int volume_points, float * rho, int r_points, double size);
+
+double modelcalculations_calculateIq_2D_3Darray(float * rho, int r_points, double step, double q, double phi);
+int modelcalculations_calculatePairCorrelation_2D_3Darray(SpacePoint * points, int volume_points, float * rho, int r_points, double size);
+SpacePoint modelcalculations_rotate(SpacePoint p, double theta, double phi, double omega);
+
+int modelcalculations_calculatePairCorrelation_2D_vector(SpacePoint * points, int volume_points, float * rho, 
+			int r_points, double bin_width, double theta_beam, double phi_beam, double omega_beam);
+
+#endif
Index: sansmodels/src/sans/models/prototypes/src/simcylinder.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/simcylinder.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/simcylinder.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,334 @@
+#include "simcylinder.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include <memory.h>
+#include "modelCalculations.h"
+#include "libCylinder.h"
+
+double test_analytical_2D(SimCylinderParameters *pars, double q, double phi) {
+	double cyl_x, cyl_y, cyl_z;
+	double q_x, q_y, q_z, lenq;
+	double theta, alpha, f, vol, sin_val, cos_val;
+	double answer;
+        
+    // Cylinder orientation
+    cyl_x = sin(pars->theta) * cos(pars->phi);
+    cyl_y = sin(pars->theta) * sin(pars->phi);
+    cyl_z = cos(pars->theta);
+     
+    // q orientation
+    
+    q_x = cos(phi);
+    q_y = sin(phi);
+    q_z = 0;
+    
+    // Length of q vector
+    //lenq = sqrt(q_x*q_x + q_y*q_y + q_z*q_z);
+    
+    // Normalize unit vector q
+    //q_x = q_x/lenq;
+    //q_y = q_y/lenq;
+    //q_z = q_z/lenq;
+    
+    // Compute the angle btw vector q and the
+    // axis of the cylinder
+    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+    
+    // The following test should always pass
+    if (fabs(cos_val)>1.0) {
+    	printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+     	return 0;
+    }
+    
+	alpha = acos( cos_val );
+	
+	// Call the IGOR library function to get the kernel
+	answer = acos(-1.0)/2.0 * CylKernel(q, pars->radius, pars->length/2.0, alpha)/sin(alpha);
+	
+	//normalize by cylinder volume
+	//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
+    vol = acos(-1.0) * pars->radius * pars->radius * pars->length;
+	answer *= vol;
+	
+	//Scale
+	answer *= pars->scale;
+	
+	//convert to [cm-1]
+	//answer *= 1.0e8;
+	
+	return answer;
+}
+
+
+/// 1D scattering function
+double simcylinder_analytical_1D(SimCylinderParameters *pars, double q) {
+	/***
+	 * Things to keep here:
+	 * 	- volume calc
+	 *  - point generation
+	 *  
+	 */
+	
+	// Check if Rho array is available
+	int volume_points;
+	int r_points;
+	int ptsGenerated;
+	
+	double r_step;
+	double vol;
+	double bin_width;
+	double twopiq;
+	double tmp;
+		
+	//printf("returning %g %g\n",test_analytical_2D(pars, q, 0.0), pars->length);
+	return test_analytical_2D(pars, q, 0.0);
+	
+	
+	// These should be parameters
+	vol = 1.0*acos(-1.0)*pars->radius*pars->radius*pars->length;
+	
+	
+	r_points      = pars->calcPars.r_points;
+	volume_points = pars->calcPars.volume_points;
+	
+	if(pars->calcPars.isPointMemAllocated==0) {
+		
+		// Call modelCalc function here to init_volume
+		
+		//twopiq        = (2*acos(-1.0)*pars->qmax);
+		//tmp           = twopiq*twopiq*twopiq*vol;
+		//volume_points = (int) floor(tmp);
+		//r_points = (int) floor(10.0*pow(tmp,0.3333));
+		//r_points      = (int) floor(pow(tmp,0.3333));
+		//printf("v, r = %d, %d\n",volume_points, r_points); 
+		
+		// TEST
+		volume_points = 1000;
+		r_points = 1000;
+		
+		pars->calcPars.volume_points = volume_points;
+		pars->calcPars.r_points      = r_points;
+		
+		// Memory allocation
+		pars->calcPars.points = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points==NULL) {
+			printf("Problem allocating memory for 1D volume points\n");
+			return -1.0;
+		}
+		pars->calcPars.isPointMemAllocated=1;
+	}
+
+    r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points;
+	
+	//printf("step=%g  r=%g  l=%g\n", r_step, pars->radius, pars->length);
+	if(pars->calcPars.isRhoAvailable==0) {
+		// Generate random points accross the volume
+		ptsGenerated = simcylinder_generatePoints(pars->calcPars.points, volume_points, pars->radius, pars->length);
+		
+		// Consistency check
+		if(ptsGenerated <= 0) {
+			// Set error code here
+			return 0;
+		}
+		
+		// Allocate memory
+		
+		pars->calcPars.rho = (double*) malloc(r_points*sizeof(double));
+		if(pars->calcPars.rho==NULL){
+			printf("Problem allocating memory for 1D correlation points\n");
+			return 0;
+		}
+		
+				
+		if(modelcalculations_calculatePairCorrelation_1D(pars->calcPars.points, volume_points, pars->calcPars.rho, r_points, r_step)<=0){
+			printf("Error occured!\n");
+			return 0;
+		};
+		pars->calcPars.isRhoAvailable=1;
+	} 
+	
+	// Calculate I(q,phi) and return that value
+	
+	return acos(-1.0)/2.0*modelcalculations_calculateIq_1D(pars->calcPars.rho, r_points, r_step, q) * vol;
+
+}
+
+/// 1D scattering function
+double simcylinder_analytical_2D(SimCylinderParameters *pars, double q, double phi) {
+	// Check if Rho array is available
+	int volume_points;
+	int r_points;
+	int ptsGenerated;
+	double bin_width;
+	double r_step;
+	double vol;
+	double retval;
+	
+	// These should be parameters
+	//r_points      = pars->calcPars.r_points;
+	//volume_points = pars->calcPars.volume_points;
+	volume_points = 50000;
+	r_points = 100;
+	if(pars->calcPars.isPointMemAllocated_2D==0) {
+		//volume_points = 1000;
+		//r_points = 100;
+		pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points_2D==NULL) {
+			printf("Problem allocating memory for 2D volume points\n");
+			return -1.0;
+		}
+		pars->calcPars.isPointMemAllocated_2D=1;
+	}
+		  
+    r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points;
+	if(pars->calcPars.isRhoAvailable_2D==0) {
+	    // Initialize random number generator 
+		int seed;    
+		seed = 10000;
+		srand(seed);
+		
+		// Generate random points accross the volume
+		ptsGenerated = simcylinder_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius, pars->length);
+		
+		// Calculate correlation function
+		pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*sizeof(float));
+		if(pars->calcPars.rho_2D==NULL){
+			printf("Problem allocating memory for 2D correlations points\n");
+			return -1.0;
+		}
+		//if(modelcalculations_calculatePairCorrelation_2D(pars->calcPars.points_2D, volume_points, pars->calcPars.rho_2D, r_points, r_step)==NULL){
+		if(modelcalculations_calculatePairCorrelation_2D_vector(pars->calcPars.points_2D, 
+			volume_points, pars->calcPars.rho_2D, r_points, r_step,
+			pars->theta, pars->phi,0.0)==0){
+			return 0;
+		};
+		pars->calcPars.isRhoAvailable_2D=1;
+	} 
+	// Calculate I(q,phi) and return that value
+	vol = acos(-1.0)*pars->radius*pars->radius*pars->length;
+	
+	//printf("in ana_2D %f %f\n",q, phi);
+	retval = modelcalculations_calculateIq_2D(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol;	
+	//printf("I=%g %f %f\n",retval, q, pars->theta);
+	return acos(-1.0)/2.0*retval;
+}
+
+/// 1D scattering function
+double simcylinder_analytical_2D_3Darray(SimCylinderParameters *pars, double q, double phi) {
+	// Check if Rho array is available
+	int volume_points;
+	int r_points;
+	int ptsGenerated;
+	double bin_width;
+	double r_step;
+	double vol;
+	
+	// These should be parameters
+	volume_points = 5000;
+	r_points = 100;
+	
+	if(pars->calcPars.isPointMemAllocated_2D==0) {
+		pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points_2D==NULL) {
+			printf("Problem allocating memory for 2D volume points\n");
+			return -1.0;
+		}
+		pars->calcPars.isPointMemAllocated_2D=1;
+	}
+		  
+    r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points;
+	
+	if(pars->calcPars.isRhoAvailable_2D==0) {
+	    // Initialize random number generator 
+		int seed;    
+		seed = 10000;
+		srand(seed);
+		 
+		// Generate random points accross the volume
+		ptsGenerated = simcylinder_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius, pars->length);
+		
+		// Calculate correlation function
+		// TODO: Check memory leaks
+		printf("A\n");
+		pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*r_points*sizeof(float));
+		if(pars->calcPars.rho_2D==NULL){
+			printf("Problem allocating memory for 2D correlations points\n");
+			return -1.0;
+		}
+		printf("B\n");
+		if(modelcalculations_calculatePairCorrelation_2D_3Darray(pars->calcPars.points_2D, 
+			volume_points, pars->calcPars.rho_2D, r_points, r_step)<0){
+			return 0;
+		};
+		printf("C\n");
+		pars->calcPars.isRhoAvailable_2D=1;
+	} 
+	// Calculate I(q,phi) and return that value
+	vol = acos(-1.0)*pars->radius*pars->radius*pars->length;
+	
+	printf("in ana_2D %f %f\n",q, phi);
+	return modelcalculations_calculateIq_2D_3Darray(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol;
+	
+	
+}
+
+
+/**
+ * Generate points randomly accross the volume
+ * @param points [SpacePoint*] Array of 3D points to be filled
+ * @param n [int] Number of points to generat
+ * @param radius [double] Radius of the sphere
+ * @return Number of points generated
+ */
+int simcylinder_generatePoints(SpacePoint * points, int n, double radius, double length) {
+	int i;
+	int testcounter;
+	double x, y, z;
+	
+	// Create points	
+	// To have a uniform density, you want to generate
+	// random points in a box and keep only those that are
+	// within the volume.
+	
+	// Initialize random number generator
+	int seed;    
+	time_t now;
+	
+	time(&now);
+	//seed = 10000;
+	
+	seed = (int)floor(fmod(now,10000));
+	//seed = 10009;
+	srand(seed);	
+	printf("Seed = %i\n", seed);
+	
+	testcounter = 0;
+		
+	memset(points,0,n*sizeof(SpacePoint));
+	for(i=0;i<n;i++) {
+		// Generate in a box centered around zero
+		x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * length/2.0;
+		
+		// reject those that are not within the volume
+		if( sqrt(x*x+y*y) < radius && fabs(z)<length/2.0) {
+			points[i].x =  x;
+			points[i].y =  y;
+			points[i].z =  z;
+			testcounter++;
+		} else {
+			i--;
+		}
+	}
+	//printf("test counter = %d\n", testcounter);
+	
+	// Consistency check
+	if(testcounter != n) {
+		return -1;
+	}
+		
+	return testcounter;
+}
Index: sansmodels/src/sans/models/prototypes/src/simcylinder.h
===================================================================
--- sansmodels/src/sans/models/prototypes/src/simcylinder.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/simcylinder.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,46 @@
+#if !defined(simcylinder_h)
+#define simcylinder_h
+
+#include "modelCalculations.h"
+
+/** Structure definition for simulated cylinder parameters 
+ * [PYTHONCLASS] = SimCylinder
+ * */
+typedef struct {
+    /// Scale factor 
+    //  [DEFAULT]=scale=1.0
+    double scale;
+
+    /// Radius of the cylinder
+    //  [DEFAULT]=radius=20.0
+    double radius;
+
+    /// Length of the cylinder
+    //  [DEFAULT]=length=400.0
+    double length;
+
+    /// Theta angle relative to beam
+    //  [DEFAULT]=theta=0.0
+    double theta;
+
+    /// Phi angle relative to x-axis
+    //  [DEFAULT]=phi=0.0
+    double phi;
+
+    /// Maximum Q (for numerical calculations)
+    //  [DEFAULT]=qmax=0.1
+    double qmax;
+    // Variables for numerical calculations 
+    CalcParameters calcPars;
+    
+} SimCylinderParameters;
+
+/// 1D scattering function
+double simcylinder_analytical_1D(SimCylinderParameters *pars, double q);
+
+/// 2D scattering function
+double simcylinder_analytical_2D(SimCylinderParameters *pars, double q, double phi);
+
+int simcylinder_generatePoints(SpacePoint * points, int n, double radius, double length);
+
+#endif
Index: sansmodels/src/sans/models/prototypes/src/simcylinder_fast.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/simcylinder_fast.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/simcylinder_fast.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,316 @@
+#include "simcylinder_fast.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include <memory.h>
+
+
+double simcylinder_fast_analytical_1D_test(SimCylinderFParameters *pars, double q) {
+	/***
+	 * Things to keep here:
+	 * 	- volume calc
+	 *  - point generation
+	 *  
+	 */
+		// Check if Rho array is available
+	int volume_points;
+	int ptsGenerated, npts;
+	double vol;
+	
+	int i, j, k, m;
+	SpacePoint *p1;
+	
+	double sum;
+	double theta_cyl, phi;
+	double phase;
+	
+	double cos_term;
+	double sin_term;
+	double qx, qy;
+	double pi_step;
+	
+	
+	volume_points = (int)floor(pars->npoints);
+	
+		
+	// Generate random points accross the volume
+	if(pars->calcPars.isPointMemAllocated_2D==0) {
+		pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points_2D==NULL) {
+			printf("Problem allocating memory for 2D volume points\n");
+			return -1.0;
+		}
+		
+		ptsGenerated = simcylinder_fast_generatePoints(pars->calcPars.points_2D, 
+			volume_points, pars->radius, pars->length, (int)floor(pars->seed));
+		
+		pars->calcPars.isPointMemAllocated_2D=1;
+	}
+		  
+	// Loop over theta_cyl
+	
+			
+	sum = 0;
+	npts = 21;
+	// Allocate temporary memory
+	p1 = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+	
+	pi_step =  acos(-1.0)/npts;
+	for(i=0; i<npts; i++) {
+		theta_cyl = pi_step * i;
+		// Rotate the points
+		
+		for(k=0; k<volume_points; k++) { 
+			p1[k] = fast_rotate(pars->calcPars.points_2D[k], theta_cyl, 0.0, 0.0);
+		}
+		
+		// Loop over phi_cyl
+		// Equivalent to looping over phi (detector)		
+		for(j=0; j<npts; j++) {
+			phi = 2*pi_step * j;
+			
+			cos_term = 0.0;
+			sin_term = 0.0;
+			
+			qx = q*cos(phi);
+			qy = q*sin(phi);
+		
+			for(m=0;m<volume_points;m++) {
+				phase = qx*p1[m].x + qy*p1[m].y;
+				cos_term += cos(phase);
+				sin_term += sin(phase);
+			}
+		
+		    sum += sin(theta_cyl)* (cos_term*cos_term + sin_term*sin_term);
+		}
+		
+	}
+	free(p1);
+	// Calculate I(q,phi) and return that value
+	vol = acos(-1.0)*pars->radius*pars->radius*pars->length;
+	sum = 1.0e8/volume_points/volume_points*vol*sum*acos(-1.0)/2;
+	return sum/npts/npts;
+}
+
+/// 1D scattering function
+double simcylinder_fast_analytical_1D(SimCylinderFParameters *pars, double q) {
+	int i, j, npts, volume_points;
+	double sum, vol, psi;
+	
+	//return simcylinder_fast_analytical_1D_test(pars, q);
+	
+	
+	sum = 0;
+	npts = 51;
+	volume_points = (int)floor(pars->npoints);
+	
+	for(i=0; i<npts; i++) {
+		psi = 2*acos(-1.0)/npts * i;
+		pars->phi = psi;
+		
+		for(j=0; j<npts; j++) {
+			pars->theta = acos(-1.0)/npts * j;
+			sum += sin(pars->theta)*simcylinder_simple_analytical_2D(pars, q, 0.0);
+		}
+	}
+	// Calculate I(q,phi) and return that value
+	vol = acos(-1.0)*pars->radius*pars->radius*pars->length;
+	sum = 1.0e8/volume_points/volume_points*vol*sum*acos(-1.0)/2;
+	return sum/npts/npts;
+}
+
+
+/// 2D scattering function
+double simcylinder_fast_analytical_2D(SimCylinderFParameters *pars, double q, double phi) {
+	int volume_points;
+	int ptsGenerated;
+	double vol;
+	
+	int i;
+	SpacePoint p1;
+	
+	double cos_term;
+	double sin_term;
+	double qx, qy, qz;
+	double phase;
+
+	volume_points = (int)floor(pars->npoints);
+	cos_term = 0.0;
+	sin_term = 0.0;
+	
+	qx = q*cos(phi);
+	qy = q*sin(phi);
+	qz = 0.0;
+		
+	// Generate random points accross the volume
+	if(pars->calcPars.isPointMemAllocated_2D==0) {
+		pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points_2D==NULL) {
+			printf("Problem allocating memory for 2D volume points\n");
+			return -1.0;
+		}
+		
+	    ptsGenerated = simcylinder_fast_generatePoints(pars->calcPars.points_2D, 
+			volume_points, pars->radius, pars->length, (int)floor(pars->seed));
+		
+		for(i=0;i<volume_points;i++) {
+			p1 = fast_rotate(pars->calcPars.points_2D[i], pars->theta, pars->phi, 0.0);
+			pars->calcPars.points_2D[i] = p1;
+		}
+				
+		pars->calcPars.isPointMemAllocated_2D=1;
+	}
+		  
+			
+	for(i=0;i<volume_points-1;i++) {
+		p1 = pars->calcPars.points_2D[i];
+		phase = qx*p1.x + qy*p1.y;
+		cos_term += cos(phase);
+		sin_term += sin(phase);
+		    
+	}
+	
+	// Calculate I(q,phi) and return that value
+	vol = acos(-1.0)*pars->radius*pars->radius*pars->length;
+	return 1.0e8/volume_points/volume_points*vol*acos(-1)/2*(cos_term*cos_term + sin_term*sin_term);	
+}
+
+double simcylinder_simple_analytical_2D(SimCylinderFParameters *pars, double q, double phi) {
+	int ptsGenerated;
+	
+	int i;
+	SpacePoint p1;
+	
+	double cos_term;
+	double sin_term;
+	double qx, qy, qz;
+	double phase;
+	int volume_points;
+	
+	cos_term = 0.0;
+	sin_term = 0.0;
+	
+	qx = q*cos(phi);
+	qy = q*sin(phi);
+	qz = 0.0;
+
+	volume_points = (int)floor(pars->npoints);
+	
+	// Generate random points accross the volume
+	if(pars->calcPars.isPointMemAllocated_2D==0) {
+		pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points_2D==NULL) {
+			printf("Problem allocating memory for 2D volume points\n");
+			return -1.0;
+		}
+		
+	    ptsGenerated = simcylinder_fast_generatePoints(pars->calcPars.points_2D, 
+			volume_points, pars->radius, pars->length, (int)floor(pars->seed));
+		
+		for(i=0;i<volume_points;i++) {
+			p1 = fast_rotate(pars->calcPars.points_2D[i], pars->theta, pars->phi, 0.0);
+			pars->calcPars.points_2D[i] = p1;
+		}
+		
+		pars->calcPars.isPointMemAllocated_2D=1;
+	}
+			
+	for(i=0;i<volume_points;i++) {
+		p1 = fast_rotate(pars->calcPars.points_2D[i], pars->theta, pars->phi, 0.0);
+		phase = qx*p1.x + qy*p1.y;
+		cos_term += cos(phase);
+		sin_term += sin(phase);
+	}
+
+    return (cos_term*cos_term + sin_term*sin_term);
+}
+
+/**
+ *  Rotation of pair correlation function
+ */
+SpacePoint fast_rotate(SpacePoint p, double theta, double phi, double omega) {
+	SpacePoint new_point;
+	double x_1, x_2;
+	double y_1, y_2;
+	double z_1, z_2;
+	
+	
+	// Omega, around z-axis (doesn't change anything for cylindrical symmetry
+	x_1 = p.x*cos(omega) - p.y*sin(omega);
+	y_1 = p.x*sin(omega) + p.y*cos(omega);
+	z_1 = p.z;	
+	
+	// Theta, around y-axis
+	x_2 = x_1*cos(theta) + z_1*sin(theta);
+	y_2 = y_1;
+	z_2 = -x_1*sin(theta) + z_1*cos(theta);
+	
+	// Phi, around z-axis
+	new_point.x = x_2*cos(phi) - y_2*sin(phi);
+	new_point.y = x_2*sin(phi) + y_2*cos(phi);
+	new_point.z = z_2;
+	
+	return new_point;
+	
+	
+}
+ 
+
+
+/**
+ * Generate points randomly accross the volume
+ * @param points [SpacePoint*] Array of 3D points to be filled
+ * @param n [int] Number of points to generat
+ * @param radius [double] Radius of the sphere
+ * @return Number of points generated
+ */
+int simcylinder_fast_generatePoints(SpacePoint * points, int n, double radius, double length, int seed) {
+	int i;
+	int testcounter;
+	double x, y, z;
+	
+	// Create points	
+	// To have a uniform density, you want to generate
+	// random points in a box and keep only those that are
+	// within the volume.
+	
+	// Initialize random number generator
+	//int seed;    
+	time_t now;
+	
+	time(&now);
+	//seed = 10000;
+	
+	//seed = (int)floor(fmod(now,10000));
+	//seed = 10009; 
+	srand(seed);	
+	//printf("Seed = %i\n", seed);
+	
+	testcounter = 0;
+		
+	memset(points,0,n*sizeof(SpacePoint));
+	for(i=0;i<n;i++) {
+		// Generate in a box centered around zero
+		x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * length/2.0;
+		
+		// reject those that are not within the volume
+		if( sqrt(x*x+y*y) < radius && fabs(z)<length/2.0) {
+			points[i].x =  x;
+			points[i].y =  y;
+			points[i].z =  z;
+			testcounter++;
+		} else {
+			i--;
+		}
+	}
+	
+	// Consistency check
+	if(testcounter != n) {
+		return -1;
+	}
+		
+	return testcounter;
+}
Index: sansmodels/src/sans/models/prototypes/src/simcylinder_fast.h
===================================================================
--- sansmodels/src/sans/models/prototypes/src/simcylinder_fast.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/simcylinder_fast.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,53 @@
+#if !defined(simcylinderf_h)
+#define simcylinderf_h
+
+#include "modelCalculations.h"
+
+/** Structure definition for simulated cylinder parameters 
+ * [PYTHONCLASS] = SimCylinderF
+ * */
+typedef struct {
+    /// Scale factor 
+    //  [DEFAULT]=scale=1.0
+    double scale;
+
+    /// Radius of the cylinder
+    //  [DEFAULT]=radius=20.0
+    double radius;
+
+    /// Length of the cylinder
+    //  [DEFAULT]=length=400.0
+    double length;
+
+    /// Theta angle relative to beam
+    //  [DEFAULT]=theta=0.0
+    double theta;
+
+    /// Phi angle relative to x-axis
+    //  [DEFAULT]=phi=0.0
+    double phi;
+
+    /// Number of space points
+    //  [DEFAULT]=npoints=50000
+    double npoints;
+
+    /// Random number seed
+    //  [DEFAULT]=seed=1000.0
+    double seed;
+
+    // Variables for numerical calculations 
+    CalcParameters calcPars;
+    
+} SimCylinderFParameters;
+
+/// 1D scattering function
+double simcylinder_fast_analytical_1D(SimCylinderFParameters *pars, double q);
+
+/// 2D scattering function
+double simcylinder_fast_analytical_2D(SimCylinderFParameters *pars, double q, double phi);
+
+int simcylinder_fast_generatePoints(SpacePoint * points, int n, double radius, double length, int seed);
+SpacePoint fast_rotate(SpacePoint p, double theta, double phi, double omega);
+double simcylinder_simple_analytical_2D(SimCylinderFParameters *pars, double q, double phi);
+
+#endif
Index: sansmodels/src/sans/models/prototypes/src/smeared_cylinder.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/smeared_cylinder.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/smeared_cylinder.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,178 @@
+/**
+ * Scattering model for a cylinder
+ * @author: Mathieu Doucet / UTK
+ */
+
+#include "cylinder.h"
+#include "smeared_cylinder.h"
+#include <math.h>
+#include "libCylinder.h"
+
+
+/**
+ * Function to evaluate 1D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+double smeared_cylinder_analytical_1D(SmearCylinderParameters *pars, double q) {
+	double dp[5];
+	int i_r;
+	double r_0, r, step_r, min_r;
+	int npts;
+	double weight, func, norm;
+	
+	// Fill paramater array
+	dp[0] = pars->scale;
+	dp[1] = pars->radius;
+	dp[2] = pars->length;
+	dp[3] = pars->contrast;
+	dp[4] = pars->background;
+	
+	if(pars->sigma_radius==0) {
+		return CylinderForm(dp, q);
+	}
+	
+	// Central value is the current value 
+	r_0     = pars->radius;
+	
+	npts = 100;
+	step_r     = 4.0*pars->sigma_radius/npts;
+	min_r     = r_0 - 2.0*pars->sigma_radius;
+	
+	
+	norm = 0.0;
+	func = 0.0;
+	
+	for (i_r=0; i_r<100; i_r++) {
+			r = min_r + step_r*i_r;
+			dp[1] = r;
+			
+			// Weigth for that position
+			weight = smeared_cylinder_dist( r, r_0, pars->sigma_radius );
+			norm += weight;
+			
+			// Evaluate I(q) at that r-value
+			func += weight * CylinderForm(dp, q);	
+	}
+	
+	return func/norm;
+}
+	
+	
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+double smeared_cylinder_analytical_2D(SmearCylinderParameters *pars, double q, double phi) {
+	CylinderParameters cyl_pars;
+	int i_theta, i_phi, i_r;
+	int n_theta, n_phi, n_r;
+	double theta_0, phi_0, r_0;
+	int npts;
+	double weight_theta, weight_phi, weight_r;
+	double min_theta, min_phi, min_r;
+	double step_theta, step_phi, step_r;
+	double func, norm;
+	double n_width = 3.0;
+	
+	// Fill paramater struct
+	cyl_pars.scale      = pars->scale;
+	cyl_pars.radius     = pars->radius;
+	cyl_pars.length     = pars->length;
+	cyl_pars.contrast   = pars->contrast;
+	cyl_pars.background = pars->background;
+	cyl_pars.cyl_phi    = pars->cyl_phi;
+	cyl_pars.cyl_theta  = pars->cyl_theta;
+	
+	theta_0 = pars->cyl_theta;
+	phi_0   = pars->cyl_phi;
+	r_0     = pars->radius;
+	
+	npts = 25;
+	step_theta = 2.0*n_width*pars->sigma_theta/npts;
+	step_phi   = 2.0*n_width*pars->sigma_phi/npts;
+	step_r     = 2.0*n_width*pars->sigma_radius/npts;
+	
+	if (step_theta>0) {
+		n_theta = npts;
+	} else {
+		n_theta = 1;
+	}
+	
+	if (step_phi>0) {
+		n_phi = npts;
+	} else {
+		n_phi = 1;
+	}
+	
+	if (step_r>0) {
+		n_r = npts;
+	} else {
+		n_r = 1;
+	}
+	
+	
+	
+	min_theta = theta_0 - n_width*pars->sigma_theta;
+	min_phi   = phi_0 - n_width*pars->sigma_phi;
+	min_r     = r_0 - n_width*pars->sigma_radius;
+	
+	func = 0.0;
+	norm = 0.0;
+	
+	for (i_theta=0; i_theta<n_theta; i_theta++) {
+			
+		// Weight for that position
+		if(pars->sigma_theta>0) {
+			cyl_pars.cyl_theta = min_theta + step_theta*i_theta;
+			weight_theta = smeared_cylinder_dist( cyl_pars.cyl_theta, theta_0, pars->sigma_theta );
+		} else {
+			weight_theta = 1.0;
+		}
+						
+		for (i_phi=0; i_phi<n_phi; i_phi++) {
+				
+			// Weight for that position
+			if(pars->sigma_phi>0) {
+				cyl_pars.cyl_phi = min_phi + step_phi*i_phi;
+				weight_phi = smeared_cylinder_dist( cyl_pars.cyl_phi, phi_0, pars->sigma_phi );
+			} else {
+				weight_phi = 1.0;
+			}
+		
+			for (i_r=0; i_r<n_r; i_r++) {
+				
+				if(pars->sigma_radius>0) {
+					cyl_pars.radius = min_r + step_r*i_r;
+					
+					// Weight for that position
+					weight_r = smeared_cylinder_dist( cyl_pars.radius, r_0, pars->sigma_radius );
+				} else {
+					weight_r = 1.0;
+				}
+					
+				// Evaluate I(q) at that r-value
+				func += weight_theta * weight_r * weight_phi * cylinder_analytical_2D(&cyl_pars, q, phi);	
+				norm += weight_theta * weight_r * weight_phi;
+			}
+		}
+	}
+	
+
+	return func/norm;
+}
+    
+double smeared_cylinder_dist( double x, double mean, double sigma ) {
+	double vary;
+	double expo;
+		
+	//return 1.0;
+		
+	vary = x-mean;
+    expo = -vary*vary/(2.0*sigma*sigma);
+	return exp(expo);
+
+}
Index: sansmodels/src/sans/models/prototypes/src/smeared_cylinder.h
===================================================================
--- sansmodels/src/sans/models/prototypes/src/smeared_cylinder.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/smeared_cylinder.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,56 @@
+#if !defined(smeared_cylinder_h)
+#define smeared_cylinder_h
+
+/** Structure definition for cylinder parameters 
+ * [PYTHONCLASS] = SmearCylinderModel
+ * */
+typedef struct {
+    /// Scale factor 
+    //  [DEFAULT]=scale=1.0
+    double scale;
+    /// Radius of the cylinder [A]
+    //  [DEFAULT]=radius=20.0 A
+    double radius;
+    /// Length of the cylinder [A]
+    //  [DEFAULT]=length=400.0 A
+    double length;
+    /// Contrast [A-2]
+    //  [DEFAULT]=contrast=3.0e-6 A-2
+    double contrast;
+	/// Incoherent Background (cm-1) 0.000
+	//  [DEFAULT]=background=0 cm-1
+	double background;    
+    /// Orientation of the cylinder axis w/respect incoming beam [rad]
+    //  [DEFAULT]=cyl_theta=1.0 rad
+    double cyl_theta;
+    /// Orientation of the cylinder in the plane of the detector [rad]
+    //  [DEFAULT]=cyl_phi=1.0 rad
+    
+    ///
+    //  [DEFAULT]=cyl_phi=0.0 rad
+    double cyl_phi;    
+    
+    ///
+    //  [DEFAULT]=sigma_theta=0.0 rad
+    double sigma_theta;
+    
+    ///
+    //  [DEFAULT]=sigma_phi=0.0 rad
+    double sigma_phi;
+    
+    ///
+    //  [DEFAULT]=sigma_radius=0.0 rad
+    double sigma_radius;
+    
+} SmearCylinderParameters;
+
+
+
+/// 1D scattering function
+double smeared_cylinder_analytical_1D(SmearCylinderParameters *pars, double q);
+
+/// 2D scattering function
+double smeared_cylinder_analytical_2D(SmearCylinderParameters *pars, double q, double phi);
+double smeared_cylinder_dist( double x, double mean, double sigma );
+
+#endif
Index: sansmodels/src/sans/models/prototypes/src/sphere_fast.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/sphere_fast.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/sphere_fast.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,160 @@
+#include "sphere_fast.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include <memory.h>
+
+
+/// 1D scattering function
+double sphere_fast_analytical_1D(SimSphereFParameters *pars, double q) {
+	int j, npts, volume_points;
+	double sum, phi, vol;
+	
+	return simcylinder_fast_analytical_2D(pars, q, 0.0);
+	
+	//return simcylinder_fast_analytical_1D_average(pars, q);
+	
+	sum = 0;
+	npts = 21;
+	volume_points = (int)floor(pars->npoints);
+	
+	for(j=0; j<npts; j++) {
+		phi = acos(-1.0)/npts * j;
+		//sum += simcylinder_simple_analytical_2D(pars, q, phi, volume_points);
+		sum += simcylinder_fast_analytical_2D(pars, q, phi);
+	}
+
+	// Calculate I(q,phi) and return that value
+	vol = 4.0/3.0*acos(-1.0)*pars->radius*pars->radius*pars->radius;
+	sum = 1.0e8/volume_points/volume_points*vol*sum*acos(-1.0)/2;
+	return sum/npts;
+}
+
+
+/// 1D scattering function
+double sphere_fast_analytical_2D(SimSphereFParameters *pars, double q, double phi) {
+	// Check if Rho array is available
+	int volume_points;
+	int r_points;
+	int ptsGenerated;
+	double bin_width;
+	double r_step;
+	double vol;
+	double retval;
+	
+	int i,j;
+	int ix_bin, iy_bin;
+	SpacePoint p1, p2;
+	
+	clock_t start;
+	clock_t finish;
+	double cos_term;
+	double sin_term;
+	double qx, qy, qz;
+	double phase;
+	double cyl_x, cyl_y, cyl_z;
+	double q_x, q_y, q_z;
+	double cos_val, alpha;
+	
+	
+	volume_points = (int)floor(pars->npoints);
+	cos_term = 0.0;
+	sin_term = 0.0;
+	
+	qx = q*cos(phi);
+	qy = q*sin(phi);
+	qz = 0.0;
+		
+	// Generate random points accross the volume
+	if(pars->calcPars.isPointMemAllocated_2D==0) {
+		pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points_2D==NULL) {
+			printf("Problem allocating memory for 2D volume points\n");
+			return -1.0;
+		}
+		
+	    ptsGenerated = sphere_fast_generatePoints(pars->calcPars.points_2D, 
+			volume_points, pars->radius, (int)floor(pars->seed));
+		
+		pars->calcPars.isPointMemAllocated_2D=1;
+	}
+		  
+			
+	for(i=0;i<volume_points-1;i++) {
+		
+		//p1 = fast_rotate(pars->calcPars.points_2D[i], pars->theta, pars->phi, 0.0);
+		p1 = pars->calcPars.points_2D[i];
+		phase = qx*p1.x + qy*p1.y;
+		cos_term += cos(phase);
+		sin_term += sin(phase);
+		    
+	}
+	
+	// Calculate I(q,phi) and return that value
+	vol = 4.0/3.0*acos(-1.0)*pars->radius*pars->radius*pars->radius;
+	
+    return 1.0e8/volume_points/volume_points*vol*(cos_term*cos_term + sin_term*sin_term);	
+    //return 1.0e8/volume_points/volume_points*vol*(cos_term*cos_term);	
+}
+
+
+ 
+
+
+/**
+ * Generate points randomly accross the volume
+ * @param points [SpacePoint*] Array of 3D points to be filled
+ * @param n [int] Number of points to generat
+ * @param radius [double] Radius of the sphere
+ * @return Number of points generated
+ */
+int sphere_fast_generatePoints(SpacePoint * points, int n, double radius, int seed) {
+	int i;
+	int testcounter;
+	double x, y, z;
+	
+	// Create points	
+	// To have a uniform density, you want to generate
+	// random points in a box and keep only those that are
+	// within the volume.
+	
+	// Initialize random number generator
+	//int seed;    
+	time_t now;
+	
+	time(&now);
+	//seed = 10000;
+	
+	//seed = (int)floor(fmod(now,10000));
+	//seed = 10009; 
+	srand(seed);	
+	//printf("Seed = %i\n", seed);
+	
+	testcounter = 0;
+		
+	memset(points,0,n*sizeof(SpacePoint));
+	for(i=0;i<n;i++) {
+		// Generate in a box centered around zero
+		x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		
+		// reject those that are not within the volume
+		if( sqrt(x*x+y*y+z*z) < radius ) {
+			points[i].x =  x;
+			points[i].y =  y;
+			points[i].z =  z;
+			testcounter++;
+		} else {
+			i--;
+		}
+	}
+	
+	// Consistency check
+	if(testcounter != n) {
+		return -1;
+	}
+		
+	return testcounter;
+}
Index: sansmodels/src/sans/models/prototypes/src/sphere_fast.h
===================================================================
--- sansmodels/src/sans/models/prototypes/src/sphere_fast.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/sphere_fast.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,39 @@
+#if !defined(simspheref_h)
+#define simspheref_h
+
+#include "modelCalculations.h"
+
+/** Structure definition for simulated cylinder parameters 
+ * [PYTHONCLASS] = SimSphereF
+ * */
+typedef struct {
+    /// Scale factor 
+    //  [DEFAULT]=scale=1.0
+    double scale;
+
+    /// Radius of the cylinder
+    //  [DEFAULT]=radius=20.0
+    double radius;
+
+    /// Number of space points
+    //  [DEFAULT]=npoints=50000
+    double npoints;
+
+    /// Random number seed
+    //  [DEFAULT]=seed=1000
+    double seed;
+
+    // Variables for numerical calculations 
+    CalcParameters calcPars;
+    
+} SimSphereFParameters;
+
+/// 1D scattering function
+double sphere_fast_analytical_1D(SimSphereFParameters *pars, double q);
+
+/// 2D scattering function
+double sphere_fast_analytical_2D(SimSphereFParameters *pars, double q, double phi);
+
+int sphere_fast_generatePoints(SpacePoint * points, int n, double radius, int seed);
+
+#endif
Index: sansmodels/src/sans/models/prototypes/src/testsphere.c
===================================================================
--- sansmodels/src/sans/models/prototypes/src/testsphere.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/testsphere.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,271 @@
+#include "testsphere.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include <memory.h>
+#include "modelCalculations.h"
+
+
+/// 1D scattering function
+double testsphere_analytical_1D(TestSphereParameters *pars, double q) {
+	/***
+	 * Things to keep here:
+	 * 	- volume calc
+	 *  - point generation
+	 *  
+	 */
+	
+	// Check if Rho array is available
+	int volume_points;
+	int r_points;
+	int ptsGenerated;
+	
+	double r_step;
+	double vol;
+	double bin_width;
+	double twopiq;
+	double tmp;
+		
+	// These should be parameters
+	vol = 4.0*acos(-1.0)/3.0*pars->radius*pars->radius*pars->radius;
+	
+	
+	r_points      = pars->calcPars.r_points;
+	volume_points = pars->calcPars.volume_points;
+	
+	if(pars->calcPars.isPointMemAllocated==0) {
+		
+		// Call modelCalc function here to init_volume
+		
+		twopiq        = (2*acos(-1.0)*pars->qmax);
+		tmp           = twopiq*twopiq*twopiq*vol;
+		volume_points = (int) floor(tmp);
+		//r_points = (int) floor(10.0*pow(tmp,0.3333));
+		r_points      = (int) floor(pow(tmp,0.3333));
+		printf("v, r = %d, %d\n",volume_points, r_points); 
+		
+		// TEST
+		volume_points = 1000;
+		r_points = 1000;
+		
+		pars->calcPars.volume_points = volume_points;
+		pars->calcPars.r_points      = r_points;
+		
+		// Memory allocation
+		pars->calcPars.points = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points==NULL) {
+			printf("Problem allocating memory for 1D volume points\n");
+			return -1.0;
+		}
+		pars->calcPars.isPointMemAllocated=1;
+	}
+	
+	if(pars->calcPars.isRhoAvailable==0) {
+		// Generate random points accross the volume
+		ptsGenerated = testsphere_generatePoints(pars->calcPars.points, volume_points, pars->radius);
+		
+		// Consistency check
+		if(ptsGenerated <= 0) {
+			// Set error code here
+			return 0;
+		}
+		
+		// Allocate memory
+		
+		pars->calcPars.rho = (double*) malloc(r_points*sizeof(double));
+		if(pars->calcPars.rho==NULL){
+			printf("Problem allocating memory for 1D correlation points\n");
+			return 0;
+		}
+		
+		
+		// Lump all this in modelCalc function
+		
+		bin_width = 2.0*pars->radius/r_points;
+		
+		if(modelcalculations_calculatePairCorrelation_1D(pars->calcPars.points, volume_points, pars->calcPars.rho, r_points, bin_width)<=0){
+			printf("Error occured!\n");
+			return 0;
+		};
+		pars->calcPars.isRhoAvailable=1;
+	} 
+	
+	// Calculate I(q,phi) and return that value
+	r_step = 2.0*pars->radius/((double)(r_points));
+	
+	return modelcalculations_calculateIq_1D(pars->calcPars.rho, r_points, r_step, q) * vol;
+
+}
+/// 1D scattering function
+double testsphere_analytical_2D_3Darray(TestSphereParameters *pars, double q, double phi) {
+	// Check if Rho array is available
+	int volume_points;
+	int r_points;
+	int ptsGenerated;
+	double bin_width;
+	double r_step;
+	double vol;
+	
+	// These should be parameters
+	volume_points = 1000;
+	r_points = 10;
+	
+	if(pars->calcPars.isPointMemAllocated_2D==0) {
+		pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points_2D==NULL) {
+			printf("Problem allocating memory for 2D volume points\n");
+			return -1.0;
+		}
+		pars->calcPars.isPointMemAllocated_2D=1;
+	}
+		  
+	//r_step = 2.0*pars->radius/((double)(r_points));
+	// Allow negative values....
+	r_step = 4.0*pars->radius/((double)(r_points));
+	
+	if(pars->calcPars.isRhoAvailable_2D==0) {
+	    // Initialize random number generator 
+		int seed;    
+		seed = 10000;
+		srand(seed);
+		 
+		// Generate random points accross the volume
+		ptsGenerated = testsphere_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius);
+		
+		// Calculate correlation function
+		pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*r_points*sizeof(float));
+		if(pars->calcPars.rho_2D==NULL){
+			printf("Problem allocating memory for 2D correlations points\n");
+			return -1.0;
+		}
+		if(modelcalculations_calculatePairCorrelation_2D_3Darray(pars->calcPars.points_2D, volume_points, pars->calcPars.rho_2D, r_points, r_step)==NULL){
+			return 0;
+		};
+		pars->calcPars.isRhoAvailable_2D=1;
+	} 
+	// Calculate I(q,phi) and return that value
+	vol = 4.0*acos(-1.0)/3.0*pars->radius*pars->radius*pars->radius;
+	
+	//printf("in ana_2D %f %f\n",q, phi);
+	return modelcalculations_calculateIq_2D_3Darray(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol;
+	
+	
+}
+
+/// 1D scattering function
+double testsphere_analytical_2D(TestSphereParameters *pars, double q, double phi) {
+	// Check if Rho array is available
+	int volume_points;
+	int r_points;
+	int ptsGenerated;
+	double bin_width;
+	double r_step;
+	double vol;
+	
+	// These should be parameters
+	//r_points      = pars->calcPars.r_points;
+	//volume_points = pars->calcPars.volume_points;
+	volume_points = 5000;
+	r_points = 100;
+	
+	if(pars->calcPars.isPointMemAllocated_2D==0) {
+		//volume_points = 2000;
+		//r_points = 100;
+		//pars->calcPars.volume_points = volume_points;
+		//pars->calcPars.r_points      = r_points;
+		pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
+		if(pars->calcPars.points_2D==NULL) {
+			printf("Problem allocating memory for 2D volume points\n");
+			return -1.0;
+		}
+		pars->calcPars.isPointMemAllocated_2D=1;
+	}
+		  
+	r_step = 2.0*pars->radius/((double)(r_points));
+	
+	if(pars->calcPars.isRhoAvailable_2D==0) {
+	    // Initialize random number generator 
+		int seed;    
+		seed = 10000;
+		srand(seed);
+		
+		// Generate random points accross the volume
+		ptsGenerated = testsphere_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius);
+		
+		// Calculate correlation function
+		pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*sizeof(float));
+		if(pars->calcPars.rho_2D==NULL){
+			printf("Problem allocating memory for 2D correlations points\n");
+			return -1.0;
+		}
+		if(modelcalculations_calculatePairCorrelation_2D(pars->calcPars.points_2D, volume_points, pars->calcPars.rho_2D, r_points, r_step)==NULL){
+		//if(modelcalculations_calculatePairCorrelation_2D_vector(pars->calcPars.points_2D, volume_points, pars->calcPars.rho_2D, r_points, r_step,0.0,0.0,0.0)==NULL){
+			return 0;
+		};
+		pars->calcPars.isRhoAvailable_2D=1;
+	} 
+	// Calculate I(q,phi) and return that value
+	vol = 4.0*acos(-1.0)/3.0*pars->radius*pars->radius*pars->radius;
+	
+	//printf("in ana_2D %f %f\n",q, phi);
+	return modelcalculations_calculateIq_2D(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol;	
+	
+}
+
+/// 1D scattering function
+double testsphere_numerical_1D(TestSphereParameters *pars, int *array, double q) { return 0;}
+
+/// 2D scattering function
+double testsphere_numerical_2D(TestSphereParameters *pars, int *array, double q, double phi) {return 0;}
+
+/**
+ * Generate points randomly accross the volume
+ * @param points [SpacePoint*] Array of 3D points to be filled
+ * @param n [int] Number of points to generat
+ * @param radius [double] Radius of the sphere
+ * @return Number of points generated
+ */
+int testsphere_generatePoints(SpacePoint * points, int n, double radius) {
+	int i;
+	int testcounter;
+	double x, y, z;
+	
+	// Create points	
+	// To have a uniform density, you want to generate
+	// random points in a box and keep only those that are
+	// within the volume.
+	
+	// Initialize random number generator
+	int seed;    
+	seed = 10000;
+	srand(seed);	
+	
+	testcounter = 0;
+		
+	memset(points,0,n*sizeof(SpacePoint));
+	for(i=0;i<n;i++) {
+		// Generate in a box centered around zero
+		x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
+		
+		// reject those that are not within the volume
+		if( sqrt(x*x+y*y+z*z) <= radius ) {
+			points[i].x =  x;
+			points[i].y =  y;
+			points[i].z =  z;
+			testcounter++;
+		} else {
+			i--;
+		}
+	}
+	//printf("test counter = %d\n", testcounter);
+	
+	// Consistency check
+	if(testcounter != n) {
+		return -1;
+	}
+		
+	return testcounter;
+}
Index: sansmodels/src/sans/models/prototypes/src/testsphere.h
===================================================================
--- sansmodels/src/sans/models/prototypes/src/testsphere.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/sans/models/prototypes/src/testsphere.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
@@ -0,0 +1,38 @@
+#if !defined(testsphere_h)
+#define testsphere_h
+
+#include "modelCalculations.h"
+
+/** Structure definition for sphere parameters 
+ * [PYTHONCLASS] = TestSphere2
+ * */
+typedef struct {
+    /// Scale factor 
+    //  [DEFAULT]=scale=1.0
+    double scale;
+    /// Radius of the sphere
+    //  [DEFAULT]=radius=20.0
+    double radius;
+    /// Maximum Q (for numerical calculations)
+    //  [DEFAULT]=qmax=0.1
+    double qmax;
+    // Variables for numerical calculations 
+    CalcParameters calcPars;
+    
+} TestSphereParameters;
+
+/// 1D scattering function
+double testsphere_analytical_1D(TestSphereParameters *pars, double q);
+
+/// 2D scattering function
+double testsphere_analytical_2D(TestSphereParameters *pars, double q, double phi);
+
+/// 1D scattering function
+double testsphere_numerical_1D(TestSphereParameters *pars, int *array, double q);
+
+/// 2D scattering function
+double testsphere_numerical_2D(TestSphereParameters *pars, int *array, double q, double phi);
+
+int testsphere_generatePoints(SpacePoint * points, int n, double radius);
+
+#endif
