Index: sansmodels/src/c_extensions/sc.c
===================================================================
--- sansmodels/src/c_extensions/sc.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,139 +1,0 @@
-/*
- * Scattering model for a SC_ParaCrystal
- */
-#include "sc.h"
-#include "libSphere.h"
-#include <math.h>
-#include <stdio.h>
-
-
-/**
- * Function to evaluate 1D scattering function
- * @param pars: parameters of the SC_ParaCrystal
- * @param q: q-value
- * @return: function value
- */
-double sc_analytical_1D(SCParameters *pars, double q) {
-	double dp[7];
-	double result;
-
-	dp[0] = pars->scale;
-	dp[1] = pars->dnn;
-	dp[2] = pars->d_factor;
-	dp[3] = pars->radius;
-	dp[4] = pars->sldSph;
-	dp[5] = pars->sldSolv;
-	dp[6] = pars->background;
-
-	result = SC_ParaCrystal(dp, q);
-	// This FIXES a singualrity the kernel in libigor.
-	if ( result == INFINITY || result == NAN){
-		result = pars->background;
-	}
-	return result;
-}
-
-/**
- * Function to evaluate 2D scattering function
- * @param pars: parameters of the SC_ParaCrystal
- * @param q: q-value
- * @return: function value
- */
-double sc_analytical_2DXY(SCParameters *pars, double qx, double qy){
-	double q;
-	q = sqrt(qx*qx+qy*qy);
-	return sc_analytical_2D_scaled(pars, q, qx/q, qy/q);
-}
-
-double sc_analytical_2D(SCParameters *pars, double q, double phi) {
-	return sc_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
-}
-
-/**
- * Function to evaluate 2D scattering function
- * @param pars: parameters of the SCCrystalModel
- * @param q: q-value
- * @param q_x: q_x / q
- * @param q_y: q_y / q
- * @return: function value
- */
-double sc_analytical_2D_scaled(SCParameters *pars, double q, double q_x, double q_y) {
-	double a3_x, a3_y, a3_z, a2_x, a2_y, a1_x, a1_y;
-	double q_z;
-	double alpha, cos_val_a3, cos_val_a2, cos_val_a1;
-	double a1_dot_q, a2_dot_q,a3_dot_q;
-	double answer;
-	double Pi = 4.0*atan(1.0);
-	double aa, Da, qDa_2, latticeScale, Zq;
-
-	double dp[5];
-  //convert angle degree to radian
-  double theta = pars->theta * Pi/180.0;
-  double phi = pars->phi * Pi/180.0;
-  double psi = pars->psi * Pi/180.0;
-	dp[0] = 1.0;
-	dp[1] = pars->radius;
-	dp[2] = pars->sldSph;
-	dp[3] = pars->sldSolv;
-	dp[4] = 0.0;
-
-
-	aa = pars->dnn;
-	Da = pars->d_factor*aa;
-	qDa_2 = pow(q*Da,2.0);
-
-	latticeScale = (4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa,3.0);
-	/// Angles here are respect to detector coordinate instead of against q coordinate(PRB 36, 3, 1754)
-    // a3 axis orientation
-    a3_x = sin(theta) * cos(phi);//negative sign here???
-    a3_y = sin(theta) * sin(phi);
-    a3_z = cos(theta);
-
-    // q vector
-    q_z = 0.0;
-
-    // Compute the angle btw vector q and the a3 axis
-    cos_val_a3 = a3_x*q_x + a3_y*q_y + a3_z*q_z;
-    alpha = acos(cos_val_a3);
-    //alpha = acos(cos_val_a3);
-    a3_dot_q = aa*q*cos_val_a3;
-    // a1 axis orientation
-    a1_x = sin(psi);
-    a1_y = cos(psi);
-
-    cos_val_a1 = a1_x*q_x + a1_y*q_y;
-    a1_dot_q = aa*q*cos_val_a1*sin(alpha);
-
-    // a2 axis orientation
-    a2_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi);
-    a2_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi);
-    // a2 axis
-    cos_val_a2 =  sin(acos(cos_val_a1));//a2_x*q_x + a2_y*q_y;
-    a2_dot_q = aa*q*cos_val_a2*sin(alpha);
-
-    // The following test should always pass
-    if (fabs(cos_val_a3)>1.0) {
-    	printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
-     	return 0;
-    }
-    // Call Zq=Z1*Z2*Z3
-    Zq = (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a1_dot_q)+exp(-qDa_2));
-    Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a2_dot_q)+exp(-qDa_2));
-    Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a3_dot_q)+exp(-qDa_2));
-
-	// Use SphereForm directly from libigor
-	answer = SphereForm(dp,q)*Zq;
-
-	//consider scales
-	answer *= latticeScale * pars->scale;
-
-	// This FIXES a singualrity the kernel in libigor.
-	if ( answer == INFINITY || answer == NAN){
-		answer = 0.0;
-	}
-
-	// add background
-	answer += pars->background;
-
-	return answer;
-}
Index: sansmodels/src/c_extensions/sc.h
===================================================================
--- sansmodels/src/c_extensions/sc.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/sc.h	(revision 6e10cffb542a4045bc3c41ac60c4abf3ed220bb2)
@@ -1,81 +1,84 @@
 #if !defined(sc_h)
 #define sc_h
+#include "parameters.hh"
 
 /**
  * Structure definition for SC_ParaCrystal parameters
  */
- //[PYTHONCLASS] = SCCrystalModel
- //[DISP_PARAMS] = radius,phi, psi, theta
- //[DESCRIPTION] =<text>P(q)=(scale/Vp)*V_lattice*P(q)*Z(q)+bkg where scale is the volume
- //					 fraction of sphere,
- //				Vp = volume of the primary particle,
- //				V_lattice = volume correction for
- //					for the crystal structure,
- //				P(q)= form factor of the sphere (normalized),
- //				Z(q)= paracrystalline structure factor
- //					for a simple cubic structure.
- //				[Simple Cubic ParaCrystal Model]
- //				Parameters;
- //				scale: volume fraction of spheres
- //				bkg:background, R: radius of sphere
- //				dnn: Nearest neighbor distance
- //				d_factor: Paracrystal distortion factor
- //				radius: radius of the spheres
- //				sldSph: SLD of the sphere
- //				sldSolv: SLD of the solvent
- //
- //		</text>
- //[FIXED]=  radius.width;phi.width;psi.width; theta.width
- //[ORIENTATION_PARAMS]= <text> phi;psi; theta; phi.width;psi.width; theta.width</text>
-
-typedef struct {
-    /// Scale factor
-    //  [DEFAULT]=scale= 1.0
-    double scale;
-
-    ///	Nearest neighbor distance [A]
-    //  [DEFAULT]=dnn=220.0 [A]
-    double dnn;
-
-    ///	Paracrystal distortion factor g
-    //  [DEFAULT]=d_factor=0.06
-    double d_factor;
-
-    ///	Radius of sphere [A]
-    //  [DEFAULT]=radius=40.0 [A]
-    double radius;
-
-    ///	sldSph [1/A^(2)]
-    //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)]
-    double sldSph;
-
-    ///	sldSolv [1/A^(2)]
-    //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)]
-    double sldSolv;
-
-	/// Incoherent Background [1/cm]
-	//  [DEFAULT]=background=0 [1/cm]
-	double background;
-    /// Orientation of the a1 axis w/respect incoming beam [deg]
-    //  [DEFAULT]=theta=0.0 [deg]
-    double theta;
-    /// Orientation of the a2 in the plane of the detector [deg]
-    //  [DEFAULT]=phi=0.0 [deg]
-    double phi;
-    /// Orientation of the a3 in the plane of the detector [deg]
-    //  [DEFAULT]=psi=0.0 [deg]
-    double psi;
-
-} SCParameters;
+//[PYTHONCLASS] = SCCrystalModel
+//[DISP_PARAMS] = radius,phi, psi, theta
+//[DESCRIPTION] =<text>P(q)=(scale/Vp)*V_lattice*P(q)*Z(q)+bkg where scale is the volume
+//					 fraction of sphere,
+//				Vp = volume of the primary particle,
+//				V_lattice = volume correction for
+//					for the crystal structure,
+//				P(q)= form factor of the sphere (normalized),
+//				Z(q)= paracrystalline structure factor
+//					for a simple cubic structure.
+//				[Simple Cubic ParaCrystal Model]
+//				Parameters;
+//				scale: volume fraction of spheres
+//				bkg:background, R: radius of sphere
+//				dnn: Nearest neighbor distance
+//				d_factor: Paracrystal distortion factor
+//				radius: radius of the spheres
+//				sldSph: SLD of the sphere
+//				sldSolv: SLD of the solvent
+//
+//		</text>
+//[FIXED]=  radius.width;phi.width;psi.width; theta.width
+//[ORIENTATION_PARAMS]= <text> phi;psi; theta; phi.width;psi.width; theta.width</text>
 
 
 
-/// 1D scattering function
-double sc_analytical_1D(SCParameters *pars, double q);
+class SCCrystalModel{
+public:
+  // Model parameters
+  /// Scale factor
+  //  [DEFAULT]=scale= 1.0
+  Parameter scale;
 
-/// 2D scattering function
-double sc_analytical_2D(SCParameters *pars, double q, double phi);
-double sc_analytical_2DXY(SCParameters *pars, double qx, double qy);
-double sc_analytical_2D_scaled(SCParameters *pars, double q, double q_x, double q_y);
+  /// Nearest neighbor distance [A]
+  //  [DEFAULT]=dnn=220.0 [A]
+  Parameter dnn;
+
+  /// Paracrystal distortion factor g
+  //  [DEFAULT]=d_factor=0.06
+  Parameter d_factor;
+
+  /// Radius of sphere [A]
+  //  [DEFAULT]=radius=40.0 [A]
+  Parameter radius;
+
+  /// sldSph [1/A^(2)]
+  //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)]
+  Parameter sldSph;
+
+  /// sldSolv [1/A^(2)]
+  //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)]
+  Parameter sldSolv;
+
+  /// Incoherent Background [1/cm]
+  //  [DEFAULT]=background=0 [1/cm]
+  Parameter background;
+  /// Orientation of the a1 axis w/respect incoming beam [deg]
+  //  [DEFAULT]=theta=0.0 [deg]
+  Parameter theta;
+  /// Orientation of the a2 in the plane of the detector [deg]
+  //  [DEFAULT]=phi=0.0 [deg]
+  Parameter phi;
+  /// Orientation of the a3 in the plane of the detector [deg]
+  //  [DEFAULT]=psi=0.0 [deg]
+  Parameter psi;
+
+  // Constructor
+  SCCrystalModel();
+
+  // Operators to get I(Q)
+  double operator()(double q);
+  double operator()(double qx, double qy);
+  double calculate_ER();
+  double evaluate_rphi(double q, double phi);
+};
 
 #endif
Index: sansmodels/src/c_models/models.hh
===================================================================
--- sansmodels/src/c_models/models.hh	(revision 37805e993ee77aefcc3df740a88439050390d7ce)
+++ sansmodels/src/c_models/models.hh	(revision 6e10cffb542a4045bc3c41ac60c4abf3ed220bb2)
@@ -27,32 +27,4 @@
 using namespace std;
 
-
-
-
-
-
-class SCCrystalModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter dnn;
-	Parameter d_factor;
-	Parameter radius;
-	Parameter sldSph;
-	Parameter sldSolv;
-	Parameter background;
-	Parameter theta;
-	Parameter phi;
-	Parameter psi;
-
-	// Constructor
-	SCCrystalModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
 
 
Index: sansmodels/src/c_models/sc.cpp
===================================================================
--- sansmodels/src/c_models/sc.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/sc.cpp	(revision 6e10cffb542a4045bc3c41ac60c4abf3ed220bb2)
@@ -25,22 +25,137 @@
 #include <stdio.h>
 using namespace std;
+#include "sc.h"
 
 extern "C" {
-	#include "libSphere.h"
-	#include "sc.h"
+#include "libSphere.h"
+}
+
+// Convenience structure
+typedef struct {
+  double scale;
+  double dnn;
+  double d_factor;
+  double radius;
+  double sldSph;
+  double sldSolv;
+  double background;
+  double theta;
+  double phi;
+  double psi;
+} SCParameters;
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the SCCrystalModel
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double sc_analytical_2D_scaled(SCParameters *pars, double q, double q_x, double q_y) {
+  double a3_x, a3_y, a3_z, a2_x, a2_y, a1_x, a1_y;
+  double q_z;
+  double alpha, cos_val_a3, cos_val_a2, cos_val_a1;
+  double a1_dot_q, a2_dot_q,a3_dot_q;
+  double answer;
+  double Pi = 4.0*atan(1.0);
+  double aa, Da, qDa_2, latticeScale, Zq;
+
+  double dp[5];
+  //convert angle degree to radian
+  double theta = pars->theta * Pi/180.0;
+  double phi = pars->phi * Pi/180.0;
+  double psi = pars->psi * Pi/180.0;
+  dp[0] = 1.0;
+  dp[1] = pars->radius;
+  dp[2] = pars->sldSph;
+  dp[3] = pars->sldSolv;
+  dp[4] = 0.0;
+
+
+  aa = pars->dnn;
+  Da = pars->d_factor*aa;
+  qDa_2 = pow(q*Da,2.0);
+
+  latticeScale = (4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa,3.0);
+  /// Angles here are respect to detector coordinate instead of against q coordinate(PRB 36, 3, 1754)
+  // a3 axis orientation
+  a3_x = sin(theta) * cos(phi);//negative sign here???
+  a3_y = sin(theta) * sin(phi);
+  a3_z = cos(theta);
+
+  // q vector
+  q_z = 0.0;
+
+  // Compute the angle btw vector q and the a3 axis
+  cos_val_a3 = a3_x*q_x + a3_y*q_y + a3_z*q_z;
+  alpha = acos(cos_val_a3);
+  //alpha = acos(cos_val_a3);
+  a3_dot_q = aa*q*cos_val_a3;
+  // a1 axis orientation
+  a1_x = sin(psi);
+  a1_y = cos(psi);
+
+  cos_val_a1 = a1_x*q_x + a1_y*q_y;
+  a1_dot_q = aa*q*cos_val_a1*sin(alpha);
+
+  // a2 axis orientation
+  a2_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi);
+  a2_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi);
+  // a2 axis
+  cos_val_a2 =  sin(acos(cos_val_a1));//a2_x*q_x + a2_y*q_y;
+  a2_dot_q = aa*q*cos_val_a2*sin(alpha);
+
+  // The following test should always pass
+  if (fabs(cos_val_a3)>1.0) {
+    printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+  // Call Zq=Z1*Z2*Z3
+  Zq = (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a1_dot_q)+exp(-qDa_2));
+  Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a2_dot_q)+exp(-qDa_2));
+  Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a3_dot_q)+exp(-qDa_2));
+
+  // Use SphereForm directly from libigor
+  answer = SphereForm(dp,q)*Zq;
+
+  //consider scales
+  answer *= latticeScale * pars->scale;
+
+  // This FIXES a singualrity the kernel in libigor.
+  if ( answer == INFINITY || answer == NAN){
+    answer = 0.0;
+  }
+
+  // add background
+  answer += pars->background;
+
+  return answer;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the SC_ParaCrystal
+ * @param q: q-value
+ * @return: function value
+ */
+static double sc_analytical_2DXY(SCParameters *pars, double qx, double qy){
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return sc_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 SCCrystalModel :: SCCrystalModel() {
-	scale      = Parameter(1.0);
-	dnn		= Parameter(220.0);
-	d_factor = Parameter(0.06);
-	radius     = Parameter(40.0, true);
-	radius.set_min(0.0);
-	sldSph   = Parameter(3.0e-6);
-	sldSolv   = Parameter(6.3e-6);
-	background = Parameter(0.0);
-	theta  = Parameter(0.0, true);
-	phi    = Parameter(0.0, true);
-	psi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  dnn		= Parameter(220.0);
+  d_factor = Parameter(0.06);
+  radius     = Parameter(40.0, true);
+  radius.set_min(0.0);
+  sldSph   = Parameter(3.0e-6);
+  sldSolv   = Parameter(6.3e-6);
+  background = Parameter(0.0);
+  theta  = Parameter(0.0, true);
+  phi    = Parameter(0.0, true);
+  psi    = Parameter(0.0, true);
 }
 
@@ -52,49 +167,49 @@
  */
 double SCCrystalModel :: operator()(double q) {
-	double dp[7];
-
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = dnn();
-	dp[2] = d_factor();
-	dp[3] = radius();
-	dp[4] = sldSph();
-	dp[5] = sldSolv();
-	dp[6] = 0.0;
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	radius.get_weights(weights_rad);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-	double result;
-
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp[3] = weights_rad[i].value;
-
-		//Un-normalize SphereForm by volume
-		result = SC_ParaCrystal(dp, q) * pow(weights_rad[i].value,3);
-		// This FIXES a singualrity the kernel in libigor.
-		if ( result == INFINITY || result == NAN){
-			result = 0.0;
-		}
-		sum += weights_rad[i].weight
-			* result;
-		//Find average volume
-		vol += weights_rad[i].weight
-			* pow(weights_rad[i].value,3);
-
-		norm += weights_rad[i].weight;
-	}
-
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-	return sum/norm + background();
+  double dp[7];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = dnn();
+  dp[2] = d_factor();
+  dp[3] = radius();
+  dp[4] = sldSph();
+  dp[5] = sldSolv();
+  dp[6] = 0.0;
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+  double result;
+
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[3] = weights_rad[i].value;
+
+    //Un-normalize SphereForm by volume
+    result = SC_ParaCrystal(dp, q) * pow(weights_rad[i].value,3);
+    // This FIXES a singualrity the kernel in libigor.
+    if ( result == INFINITY || result == NAN){
+      result = 0.0;
+    }
+    sum += weights_rad[i].weight
+        * result;
+    //Find average volume
+    vol += weights_rad[i].weight
+        * pow(weights_rad[i].value,3);
+
+    norm += weights_rad[i].weight;
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+  return sum/norm + background();
 }
 
@@ -106,90 +221,90 @@
  */
 double SCCrystalModel :: operator()(double qx, double qy) {
-	SCParameters dp;
-	dp.scale      = scale();
-	dp.dnn   = dnn();
-	dp.d_factor   = d_factor();
-	dp.radius  = radius();
-	dp.sldSph   = sldSph();
-	dp.sldSolv   = sldSolv();
-	dp.background = 0.0;
-	dp.theta  = theta();
-	dp.phi    = phi();
-	dp.psi    = psi();
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	radius.get_weights(weights_rad);
-
-	// Get angular averaging for theta
-	vector<WeightPoint> weights_theta;
-	theta.get_weights(weights_theta);
-
-	// Get angular averaging for phi
-	vector<WeightPoint> weights_phi;
-	phi.get_weights(weights_phi);
-
-	// Get angular averaging for psi
-	vector<WeightPoint> weights_psi;
-	psi.get_weights(weights_psi);
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double norm_vol = 0.0;
-	double vol = 0.0;
-	double pi = 4.0*atan(1.0);
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp.radius = weights_rad[i].value;
-		// Average over theta distribution
-		for(size_t j=0; j< weights_theta.size(); j++) {
-			dp.theta = weights_theta[j].value;
-			// Average over phi distribution
-			for(size_t k=0; k< weights_phi.size(); k++) {
-				dp.phi = weights_phi[k].value;
-				// Average over phi distribution
-				for(size_t l=0; l< weights_psi.size(); l++) {
-					dp.psi = weights_psi[l].value;
-					//Un-normalize SphereForm by volume
-					double _ptvalue = weights_rad[i].weight
-										* weights_theta[j].weight
-										* weights_phi[k].weight
-										* weights_psi[l].weight
-										* sc_analytical_2DXY(&dp, qx, qy);
-										//* pow(weights_rad[i].value,3.0);
-					// Consider when there is infinte or nan.
-					if ( _ptvalue == INFINITY || _ptvalue == NAN){
-						_ptvalue = 0.0;
-					}
-					if (weights_theta.size()>1) {
-						_ptvalue *= fabs(sin(weights_theta[j].value*pi/180.0));
-					}
-					sum += _ptvalue;
-					// This model dose not need the volume of spheres correction!!!
-					//Find average volume
-					//vol += weights_rad[i].weight
-					//	* pow(weights_rad[i].value,3);
-					//Find norm for volume
-					//norm_vol += weights_rad[i].weight;
-
-					norm += weights_rad[i].weight
-							* weights_theta[j].weight
-							* weights_phi[k].weight
-							* weights_psi[l].weight;
-				}
-			}
-		}
-	}
-	// Averaging in theta needs an extra normalization
-	// factor to account for the sin(theta) term in the
-	// integration (see documentation).
-	if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-	if (vol != 0.0 && norm_vol != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm_vol);}
-
-	return sum/norm + background();
+  SCParameters dp;
+  dp.scale      = scale();
+  dp.dnn   = dnn();
+  dp.d_factor   = d_factor();
+  dp.radius  = radius();
+  dp.sldSph   = sldSph();
+  dp.sldSolv   = sldSolv();
+  dp.background = 0.0;
+  dp.theta  = theta();
+  dp.phi    = phi();
+  dp.psi    = psi();
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  phi.get_weights(weights_phi);
+
+  // Get angular averaging for psi
+  vector<WeightPoint> weights_psi;
+  psi.get_weights(weights_psi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp.radius = weights_rad[i].value;
+    // Average over theta distribution
+    for(size_t j=0; j< weights_theta.size(); j++) {
+      dp.theta = weights_theta[j].value;
+      // Average over phi distribution
+      for(size_t k=0; k< weights_phi.size(); k++) {
+        dp.phi = weights_phi[k].value;
+        // Average over phi distribution
+        for(size_t l=0; l< weights_psi.size(); l++) {
+          dp.psi = weights_psi[l].value;
+          //Un-normalize SphereForm by volume
+          double _ptvalue = weights_rad[i].weight
+              * weights_theta[j].weight
+              * weights_phi[k].weight
+              * weights_psi[l].weight
+              * sc_analytical_2DXY(&dp, qx, qy);
+          //* pow(weights_rad[i].value,3.0);
+          // Consider when there is infinte or nan.
+          if ( _ptvalue == INFINITY || _ptvalue == NAN){
+            _ptvalue = 0.0;
+          }
+          if (weights_theta.size()>1) {
+            _ptvalue *= fabs(sin(weights_theta[j].value*pi/180.0));
+          }
+          sum += _ptvalue;
+          // This model dose not need the volume of spheres correction!!!
+          //Find average volume
+          //vol += weights_rad[i].weight
+          //	* pow(weights_rad[i].value,3);
+          //Find norm for volume
+          //norm_vol += weights_rad[i].weight;
+
+          norm += weights_rad[i].weight
+              * weights_theta[j].weight
+              * weights_phi[k].weight
+              * weights_psi[l].weight;
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -202,5 +317,5 @@
  */
 double SCCrystalModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -210,5 +325,5 @@
  */
 double SCCrystalModel :: calculate_ER() {
-	//NOT implemented yet!!!
-	return 0.0;
-}
+  //NOT implemented yet!!!
+  return 0.0;
+}
Index: sansmodels/src/python_wrapper/CSCCrystalModel.cpp
===================================================================
--- sansmodels/src/python_wrapper/CSCCrystalModel.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,620 +1,0 @@
-/**
-	This software was developed by the University of Tennessee as part of the
-	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
-	project funded by the US National Science Foundation.
-
-	If you use DANSE applications to do scientific research that leads to
-	publication, we ask that you acknowledge the use of the software with the
-	following sentence:
-
-	"This work benefited from DANSE software developed under NSF award DMR-0520547."
-
-	copyright 2008, University of Tennessee
- */
-
-/** CSCCrystalModel
- *
- * C extension 
- *
- * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
- *          DO NOT MODIFY THIS FILE, MODIFY sc.h
- *          AND RE-RUN THE GENERATOR SCRIPT
- *
- */
-#define NO_IMPORT_ARRAY
-#define PY_ARRAY_UNIQUE_SYMBOL PyArray_API_sans
- 
-extern "C" {
-#include <Python.h>
-#include <arrayobject.h>
-#include "structmember.h"
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <time.h>
-#include "sc.h"
-}
-
-#include "models.hh"
-#include "dispersion_visitor.hh"
-
-/// Error object for raised exceptions
-static PyObject * CSCCrystalModelError = NULL;
-
-
-// Class definition
-typedef struct {
-    PyObject_HEAD
-    /// Parameters
-    PyObject * params;
-    /// Dispersion parameters
-    PyObject * dispersion;
-    /// Underlying model object
-    SCCrystalModel * model;
-    /// Log for unit testing
-    PyObject * log;
-} CSCCrystalModel;
-
-
-static void
-CSCCrystalModel_dealloc(CSCCrystalModel* self)
-{
-    Py_DECREF(self->params);
-    Py_DECREF(self->dispersion);
-    Py_DECREF(self->log);
-    delete self->model;
-    self->ob_type->tp_free((PyObject*)self);
-    
-
-}
-
-static PyObject *
-CSCCrystalModel_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
-{
-    CSCCrystalModel *self;
-    
-    self = (CSCCrystalModel *)type->tp_alloc(type, 0);
-   
-    return (PyObject *)self;
-}
-
-static int
-CSCCrystalModel_init(CSCCrystalModel *self, PyObject *args, PyObject *kwds)
-{
-    if (self != NULL) {
-    	
-    	// Create parameters
-        self->params = PyDict_New();
-        self->dispersion = PyDict_New();
-        self->model = new SCCrystalModel();
-        
-        // Initialize parameter dictionary
-        PyDict_SetItemString(self->params,"phi",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000000000));
-        PyDict_SetItemString(self->params,"dnn",Py_BuildValue("d",220.000000000000));
-        PyDict_SetItemString(self->params,"sldSph",Py_BuildValue("d",0.000003000000));
-        PyDict_SetItemString(self->params,"d_factor",Py_BuildValue("d",0.060000000000));
-        PyDict_SetItemString(self->params,"psi",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",40.000000000000));
-        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"theta",Py_BuildValue("d",0.000000000000));
-        PyDict_SetItemString(self->params,"sldSolv",Py_BuildValue("d",0.000006300000));
-        // Initialize dispersion / averaging parameter dict
-        DispersionVisitor* visitor = new DispersionVisitor();
-        PyObject * disp_dict;
-        disp_dict = PyDict_New();
-        self->model->radius.dispersion->accept_as_source(visitor, self->model->radius.dispersion, disp_dict);
-        PyDict_SetItemString(self->dispersion, "radius", disp_dict);
-        disp_dict = PyDict_New();
-        self->model->phi.dispersion->accept_as_source(visitor, self->model->phi.dispersion, disp_dict);
-        PyDict_SetItemString(self->dispersion, "phi", disp_dict);
-        disp_dict = PyDict_New();
-        self->model->psi.dispersion->accept_as_source(visitor, self->model->psi.dispersion, disp_dict);
-        PyDict_SetItemString(self->dispersion, "psi", disp_dict);
-        disp_dict = PyDict_New();
-        self->model->theta.dispersion->accept_as_source(visitor, self->model->theta.dispersion, disp_dict);
-        PyDict_SetItemString(self->dispersion, "theta", disp_dict);
-
-
-         
-        // Create empty log
-        self->log = PyDict_New();
-        
-        
-
-    }
-    return 0;
-}
-
-static char name_params[] = "params";
-static char def_params[] = "Parameters";
-static char name_dispersion[] = "dispersion";
-static char def_dispersion[] = "Dispersion parameters";
-static char name_log[] = "log";
-static char def_log[] = "Log";
-
-static PyMemberDef CSCCrystalModel_members[] = {
-    {name_params, T_OBJECT, offsetof(CSCCrystalModel, params), 0, def_params},
-	{name_dispersion, T_OBJECT, offsetof(CSCCrystalModel, dispersion), 0, def_dispersion},     
-    {name_log, T_OBJECT, offsetof(CSCCrystalModel, log), 0, def_log},
-    {NULL}  /* Sentinel */
-};
-
-/** Read double from PyObject
-    @param p PyObject
-    @return double
-*/
-double CSCCrystalModel_readDouble(PyObject *p) {
-    if (PyFloat_Check(p)==1) {
-        return (double)(((PyFloatObject *)(p))->ob_fval);
-    } else if (PyInt_Check(p)==1) {
-        return (double)(((PyIntObject *)(p))->ob_ival);
-    } else if (PyLong_Check(p)==1) {
-        return (double)PyLong_AsLong(p);
-    } else {
-        return 0.0;
-    }
-}
-/**
- * Function to call to evaluate model
- * @param args: input numpy array q[] 
- * @return: numpy array object 
- */
- 
-static PyObject *evaluateOneDim(SCCrystalModel* model, PyArrayObject *q){
-    PyArrayObject *result;
-   
-    // Check validity of array q , q must be of dimension 1, an array of double
-    if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE)
-    {
-        //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num;
-        //PyErr_SetString(PyExc_ValueError , message);
-        return NULL;
-    }
-    result = (PyArrayObject *)PyArray_FromDims(q->nd, (int *)(q->dimensions), PyArray_DOUBLE);
-	if (result == NULL) {
-        const char * message= "Could not create result ";
-        PyErr_SetString(PyExc_RuntimeError , message);
-		return NULL;
-	}
-#pragma omp parallel for
-	 for (int i = 0; i < q->dimensions[0]; i++){
-      double q_value  = *(double *)(q->data + i*q->strides[0]);
-      double *result_value = (double *)(result->data + i*result->strides[0]);
-      *result_value =(*model)(q_value);
-	}
-    return PyArray_Return(result); 
- }
-
- /**
- * Function to call to evaluate model
- * @param args: input numpy array  [x[],y[]]
- * @return: numpy array object 
- */
- static PyObject * evaluateTwoDimXY( SCCrystalModel* model, 
-                              PyArrayObject *x, PyArrayObject *y)
- {
-    PyArrayObject *result;
-    int x_len, y_len, dims[1];
-    //check validity of input vectors
-    if (x->nd != 1 || x->descr->type_num != PyArray_DOUBLE
-        || y->nd != 1 || y->descr->type_num != PyArray_DOUBLE
-        || y->dimensions[0] != x->dimensions[0]){
-        const char * message= "evaluateTwoDimXY  expect 2 numpy arrays";
-        PyErr_SetString(PyExc_ValueError , message); 
-        return NULL;
-    }
-   
-	if (PyArray_Check(x) && PyArray_Check(y)) {
-		
-	    x_len = dims[0]= x->dimensions[0];
-        y_len = dims[0]= y->dimensions[0];
-	    
-	    // Make a new double matrix of same dims
-        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
-        if (result == NULL){
-	    const char * message= "Could not create result ";
-        PyErr_SetString(PyExc_RuntimeError , message);
-	    return NULL;
-	    }
-       
-        /* Do the calculation. */
-#pragma omp parallel for
-        for (int i=0; i< x_len; i++) {
-            double x_value = *(double *)(x->data + i*x->strides[0]);
-  		    double y_value = *(double *)(y->data + i*y->strides[0]);
-  			double *result_value = (double *)(result->data +
-  			      i*result->strides[0]);
-  			*result_value = (*model)(x_value, y_value);
-        }           
-        return PyArray_Return(result); 
-        
-        }else{
-		    PyErr_SetString(CSCCrystalModelError, 
-                   "CSCCrystalModel.evaluateTwoDimXY couldn't run.");
-	        return NULL;
-		}      	
-}
-/**
- *  evalDistribution function evaluate a model function with input vector
- *  @param args: input q as vector or [qx, qy] where qx, qy are vectors
- *
- */ 
-static PyObject * evalDistribution(CSCCrystalModel *self, PyObject *args){
-	PyObject *qx, *qy;
-	PyArrayObject * pars;
-	int npars ,mpars;
-	
-	// Get parameters
-	
-	    // Reader parameter dictionary
-    self->model->phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "phi") );
-    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->dnn = PyFloat_AsDouble( PyDict_GetItemString(self->params, "dnn") );
-    self->model->sldSph = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSph") );
-    self->model->d_factor = PyFloat_AsDouble( PyDict_GetItemString(self->params, "d_factor") );
-    self->model->psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "psi") );
-    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
-    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
-    self->model->theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "theta") );
-    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
-    // Read in dispersion parameters
-    PyObject* disp_dict;
-    DispersionVisitor* visitor = new DispersionVisitor();
-    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
-    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "phi");
-    self->model->phi.dispersion->accept_as_destination(visitor, self->model->phi.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "psi");
-    self->model->psi.dispersion->accept_as_destination(visitor, self->model->psi.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "theta");
-    self->model->theta.dispersion->accept_as_destination(visitor, self->model->theta.dispersion, disp_dict);
-
-	
-	// Get input and determine whether we have to supply a 1D or 2D return value.
-	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
-	    PyErr_SetString(CSCCrystalModelError, 
-	    	"CSCCrystalModel.evalDistribution expects a q value.");
-		return NULL;
-	}
-    // Check params
-	
-    if(PyArray_Check(pars)==1) {
-		
-	    // Length of list should 1 or 2
-	    npars = pars->nd; 
-	    if(npars==1) {
-	        // input is a numpy array
-	        if (PyArray_Check(pars)) {
-		        return evaluateOneDim(self->model, (PyArrayObject*)pars); 
-		    }
-		}else{
-		    PyErr_SetString(CSCCrystalModelError, 
-                   "CSCCrystalModel.evalDistribution expect numpy array of one dimension.");
-	        return NULL;
-		}
-    }else if( PyList_Check(pars)==1) {
-    	// Length of list should be 2 for I(qx,qy)
-	    mpars = PyList_GET_SIZE(pars); 
-	    if(mpars!=2) {
-	    	PyErr_SetString(CSCCrystalModelError, 
-	    		"CSCCrystalModel.evalDistribution expects a list of dimension 2.");
-	    	return NULL;
-	    }
-	     qx = PyList_GET_ITEM(pars,0);
-	     qy = PyList_GET_ITEM(pars,1);
-	     if (PyArray_Check(qx) && PyArray_Check(qy)) {
-	         return evaluateTwoDimXY(self->model, (PyArrayObject*)qx,
-		           (PyArrayObject*)qy);
-		 }else{
-		    PyErr_SetString(CSCCrystalModelError, 
-                   "CSCCrystalModel.evalDistribution expect 2 numpy arrays in list.");
-	        return NULL;
-	     }
-	}
-	PyErr_SetString(CSCCrystalModelError, 
-                   "CSCCrystalModel.evalDistribution couln't be run.");
-	return NULL;
-	
-}
-
-/**
- * Function to call to evaluate model
- * @param args: input q or [q,phi]
- * @return: function value
- */
-static PyObject * run(CSCCrystalModel *self, PyObject *args) {
-	double q_value, phi_value;
-	PyObject* pars;
-	int npars;
-	
-	// Get parameters
-	
-	    // Reader parameter dictionary
-    self->model->phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "phi") );
-    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->dnn = PyFloat_AsDouble( PyDict_GetItemString(self->params, "dnn") );
-    self->model->sldSph = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSph") );
-    self->model->d_factor = PyFloat_AsDouble( PyDict_GetItemString(self->params, "d_factor") );
-    self->model->psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "psi") );
-    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
-    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
-    self->model->theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "theta") );
-    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
-    // Read in dispersion parameters
-    PyObject* disp_dict;
-    DispersionVisitor* visitor = new DispersionVisitor();
-    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
-    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "phi");
-    self->model->phi.dispersion->accept_as_destination(visitor, self->model->phi.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "psi");
-    self->model->psi.dispersion->accept_as_destination(visitor, self->model->psi.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "theta");
-    self->model->theta.dispersion->accept_as_destination(visitor, self->model->theta.dispersion, disp_dict);
-
-	
-	// Get input and determine whether we have to supply a 1D or 2D return value.
-	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
-	    PyErr_SetString(CSCCrystalModelError, 
-	    	"CSCCrystalModel.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(CSCCrystalModelError, 
-	    		"CSCCrystalModel.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 = CSCCrystalModel_readDouble(PyList_GET_ITEM(pars,0));
-	    phi_value = CSCCrystalModel_readDouble(PyList_GET_ITEM(pars,1));
-	    // Skip zero
-	    if (q_value==0) {
-	    	return Py_BuildValue("d",0.0);
-	    }
-		return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value));
-
-	} else {
-
-		// We have a scalar q, we will evaluate I(q)
-		q_value = CSCCrystalModel_readDouble(pars);		
-		
-		return Py_BuildValue("d",(*(self->model))(q_value));
-	}	
-}
-/**
- * Function to call to calculate_ER
- * @return: effective radius value 
- */
-static PyObject * calculate_ER(CSCCrystalModel *self) {
-
-	// Get parameters
-	
-	    // Reader parameter dictionary
-    self->model->phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "phi") );
-    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->dnn = PyFloat_AsDouble( PyDict_GetItemString(self->params, "dnn") );
-    self->model->sldSph = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSph") );
-    self->model->d_factor = PyFloat_AsDouble( PyDict_GetItemString(self->params, "d_factor") );
-    self->model->psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "psi") );
-    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
-    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
-    self->model->theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "theta") );
-    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
-    // Read in dispersion parameters
-    PyObject* disp_dict;
-    DispersionVisitor* visitor = new DispersionVisitor();
-    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
-    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "phi");
-    self->model->phi.dispersion->accept_as_destination(visitor, self->model->phi.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "psi");
-    self->model->psi.dispersion->accept_as_destination(visitor, self->model->psi.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "theta");
-    self->model->theta.dispersion->accept_as_destination(visitor, self->model->theta.dispersion, disp_dict);
-
-		
-	return Py_BuildValue("d",(*(self->model)).calculate_ER());
-
-}
-/**
- * Function to call to evaluate model in cartesian coordinates
- * @param args: input q or [qx, qy]]
- * @return: function value
- */
-static PyObject * runXY(CSCCrystalModel *self, PyObject *args) {
-	double qx_value, qy_value;
-	PyObject* pars;
-	int npars;
-	
-	// Get parameters
-	
-	    // Reader parameter dictionary
-    self->model->phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "phi") );
-    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") );
-    self->model->dnn = PyFloat_AsDouble( PyDict_GetItemString(self->params, "dnn") );
-    self->model->sldSph = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSph") );
-    self->model->d_factor = PyFloat_AsDouble( PyDict_GetItemString(self->params, "d_factor") );
-    self->model->psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "psi") );
-    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") );
-    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") );
-    self->model->theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "theta") );
-    self->model->sldSolv = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sldSolv") );
-    // Read in dispersion parameters
-    PyObject* disp_dict;
-    DispersionVisitor* visitor = new DispersionVisitor();
-    disp_dict = PyDict_GetItemString(self->dispersion, "radius");
-    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "phi");
-    self->model->phi.dispersion->accept_as_destination(visitor, self->model->phi.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "psi");
-    self->model->psi.dispersion->accept_as_destination(visitor, self->model->psi.dispersion, disp_dict);
-    disp_dict = PyDict_GetItemString(self->dispersion, "theta");
-    self->model->theta.dispersion->accept_as_destination(visitor, self->model->theta.dispersion, disp_dict);
-
-	
-	// Get input and determine whether we have to supply a 1D or 2D return value.
-	if ( !PyArg_ParseTuple(args,"O",&pars) ) {
-	    PyErr_SetString(CSCCrystalModelError, 
-	    	"CSCCrystalModel.run expects a q value.");
-		return NULL;
-	}
-	  
-	// Check params
-	if( PyList_Check(pars)==1) {
-		
-		// Length of list should be 2 for I(qx, qy))
-	    npars = PyList_GET_SIZE(pars); 
-	    if(npars!=2) {
-	    	PyErr_SetString(CSCCrystalModelError, 
-	    		"CSCCrystalModel.run expects a double or a list of dimension 2.");
-	    	return NULL;
-	    }
-	    // We have a vector q, get the qx and qy values at which
-	    // to evaluate I(qx,qy)
-	    qx_value = CSCCrystalModel_readDouble(PyList_GET_ITEM(pars,0));
-	    qy_value = CSCCrystalModel_readDouble(PyList_GET_ITEM(pars,1));
-	    return Py_BuildValue("d",(*(self->model))(qx_value,qy_value));
-
-	} else {
-
-		// We have a scalar q, we will evaluate I(q)
-		qx_value = CSCCrystalModel_readDouble(pars);		
-		
-		return Py_BuildValue("d",(*(self->model))(qx_value));
-	}	
-}
-
-static PyObject * reset(CSCCrystalModel *self, PyObject *args) {
-    
-
-    return Py_BuildValue("d",0.0);
-}
-
-static PyObject * set_dispersion(CSCCrystalModel *self, PyObject *args) {
-	PyObject * disp;
-	const char * par_name;
-
-	if ( !PyArg_ParseTuple(args,"sO", &par_name, &disp) ) {
-	    PyErr_SetString(CSCCrystalModelError,
-	    	"CSCCrystalModel.set_dispersion expects a DispersionModel object.");
-		return NULL;
-	}
-	void *temp = PyCObject_AsVoidPtr(disp);
-	DispersionModel * dispersion = static_cast<DispersionModel *>(temp);
-
-
-	// Ugliness necessary to go from python to C
-	    // TODO: refactor this
-    if (!strcmp(par_name, "radius")) {
-        self->model->radius.dispersion = dispersion;
-    } else    if (!strcmp(par_name, "phi")) {
-        self->model->phi.dispersion = dispersion;
-    } else    if (!strcmp(par_name, "psi")) {
-        self->model->psi.dispersion = dispersion;
-    } else    if (!strcmp(par_name, "theta")) {
-        self->model->theta.dispersion = dispersion;
-    } else {
-	    PyErr_SetString(CSCCrystalModelError,
-	    	"CSCCrystalModel.set_dispersion expects a valid parameter name.");
-		return NULL;
-	}
-
-	DispersionVisitor* visitor = new DispersionVisitor();
-	PyObject * disp_dict = PyDict_New();
-	dispersion->accept_as_source(visitor, dispersion, disp_dict);
-	PyDict_SetItemString(self->dispersion, par_name, disp_dict);
-    return Py_BuildValue("i",1);
-}
-
-
-static PyMethodDef CSCCrystalModel_methods[] = {
-    {"run",      (PyCFunction)run     , METH_VARARGS,
-      "Evaluate the model at a given Q or Q, phi"},
-    {"runXY",      (PyCFunction)runXY     , METH_VARARGS,
-      "Evaluate the model at a given Q or Qx, Qy"},
-    {"calculate_ER",      (PyCFunction)calculate_ER     , METH_VARARGS,
-      "Evaluate the model at a given Q or Q, phi"},
-      
-    {"evalDistribution",  (PyCFunction)evalDistribution , METH_VARARGS,
-      "Evaluate the model at a given Q or Qx, Qy vector "},
-    {"reset",    (PyCFunction)reset   , METH_VARARGS,
-      "Reset pair correlation"},
-    {"set_dispersion",      (PyCFunction)set_dispersion     , METH_VARARGS,
-      "Set the dispersion model for a given parameter"},
-   {NULL}
-};
-
-static PyTypeObject CSCCrystalModelType = {
-    PyObject_HEAD_INIT(NULL)
-    0,                         /*ob_size*/
-    "CSCCrystalModel",             /*tp_name*/
-    sizeof(CSCCrystalModel),             /*tp_basicsize*/
-    0,                         /*tp_itemsize*/
-    (destructor)CSCCrystalModel_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*/
-    "CSCCrystalModel objects",           /* tp_doc */
-    0,		               /* tp_traverse */
-    0,		               /* tp_clear */
-    0,		               /* tp_richcompare */
-    0,		               /* tp_weaklistoffset */
-    0,		               /* tp_iter */
-    0,		               /* tp_iternext */
-    CSCCrystalModel_methods,             /* tp_methods */
-    CSCCrystalModel_members,             /* tp_members */
-    0,                         /* tp_getset */
-    0,                         /* tp_base */
-    0,                         /* tp_dict */
-    0,                         /* tp_descr_get */
-    0,                         /* tp_descr_set */
-    0,                         /* tp_dictoffset */
-    (initproc)CSCCrystalModel_init,      /* tp_init */
-    0,                         /* tp_alloc */
-    CSCCrystalModel_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 addCSCCrystalModel(PyObject *module) {
-	PyObject *d;
-	
-    if (PyType_Ready(&CSCCrystalModelType) < 0)
-        return;
-
-    Py_INCREF(&CSCCrystalModelType);
-    PyModule_AddObject(module, "CSCCrystalModel", (PyObject *)&CSCCrystalModelType);
-    
-    d = PyModule_GetDict(module);
-    static char error_name[] = "CSCCrystalModel.error";
-    CSCCrystalModelError = PyErr_NewException(error_name, NULL, NULL);
-    PyDict_SetItemString(d, "CSCCrystalModelError", CSCCrystalModelError);
-}
-
