Index: nsmodels/src/c_extensions/onion.c
===================================================================
--- sansmodels/src/c_extensions/onion.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,267 +1,0 @@
-/**
- * Scattering model for a onion
- */
-
-#include <math.h>
-#include "onion.h"
-#include <stdio.h>
-#include <stdlib.h>
-// some details can be found in sld_cal.c
-double so_kernel(double dp[], double q) {
-	int n = dp[0];
-	double scale = dp[1];
-	double rad_core0 = dp[2];
-	double sld_core0 = dp[3];
-	double sld_solv = dp[4];
-	double background = dp[5];
-	int i,j;
-  double bes,fun,alpha,f,vol,vol_pre,vol_sub,qr,r,contr,f2;
-  double sign;
-  double pi;
-  double r0 = 0.0;
-
-  double *sld_out;
-  double *slope;
-  double *sld_in;
-  double *thick;
-  double *A;
-  int *fun_type;
-
-  sld_out = (double*)malloc((n+2)*sizeof(double));
-  slope = (double*)malloc((n+2)*sizeof(double));
-  sld_in = (double*)malloc((n+2)*sizeof(double));
-  thick = (double*)malloc((n+2)*sizeof(double));
-  A = (double*)malloc((n+2)*sizeof(double));
-  fun_type = (int*)malloc((n+2)*sizeof(int));
-
-	for (i =1; i<=n; i++){
-		sld_out[i] = dp[i+5];
-		sld_in[i] = dp[i+15];
-		A[i] = dp[i+25];
-		thick[i] = dp[i+35];
-		fun_type[i] = dp[i+45];
-	}
-	sld_out[0] = sld_core0;
-	sld_out[n+1] = sld_solv;
-	sld_in[0] = sld_core0;
-	sld_in[n+1] = sld_solv;
-	thick[0] = rad_core0;
-	thick[n+1] = 1e+10;
-	A[0] = 0.0;
-	A[n+1] = 0.0;
-	fun_type[0] = 0;
-	fun_type[n+1] = 0;
-
-
-    pi = 4.0*atan(1.0);
-    f = 0.0;
-    r = 0.0;
-    vol = 0.0;
-    vol_pre = 0.0;
-    vol_sub = 0.0;
-
-    for (i =0; i<= n+1; i++){
-    	if (thick[i] == 0.0){
-    		continue;
-    	}
-    	if (fun_type[i]== 0 ){
-    		slope[i] = 0.0;
-    		A[i] = 0.0;
-    	}
-    	vol_pre = vol;
-        switch(fun_type[i]){
-            case 2 :
-					r0 = r;
-					if (A[i] == 0.0){
-						slope[i] = 0.0;
-					}
-					else{
-						slope[i] = (sld_out[i]-sld_in[i])/(exp(A[i])-1.0);
-					}
-                    for (j=0; j<2; j++){
-                        if ( i == 0 && j == 0){
-                            continue;
-                            }
-                        if (i == n+1 && j == 1){
-                            continue;
-                            }
-                        if ( j == 1){
-                            sign = 1.0;
-                            r += thick[i];
-                            }
-                        else{
-                            sign = -1.0;
-                       		}
-                        qr = q * r;
-                        alpha = A[i] * r/thick[i];
-                        fun = 0.0;
-                        if(qr == 0.0){
-                            fun = sign * 1.0;
-                            bes = sign * 1.0;
-                            }
-                        else{
-                            if (fabs(A[i]) > 0.0 ){
-                                fun = 3.0 * ((alpha*alpha - qr * qr) * sin(qr) - 2.0 * alpha * qr * cos(qr))/ ((alpha * alpha + qr * qr) * (alpha * alpha + qr * qr) * qr);
-                                fun = fun - 3.0 * (alpha * sin(qr) - qr * cos(qr)) / ((alpha * alpha + qr * qr) * qr);
-                                fun = - sign *fun;
-                                bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
-                                }
-                            else {
-                                fun = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
-                                bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
-                                }
-                            }
-                        contr = slope[i]*exp(A[i]*(r-r0)/thick[i]);
-
-                        vol = 4.0 * pi / 3.0 * r * r * r;
-                        //if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && A[i]==0.0){
-                        //	vol_sub += (vol_pre - vol);
-                        //}
-                        f += vol * (contr * (fun) + (sld_in[i]-slope[i]) * bes);
-                        }
-                        break;
-            default :
-						if (fun_type[i]==0){
-							slope[i] = 0.0;
-						}
-						else{
-							slope[i]= (sld_out[i] -sld_in[i])/thick[i];
-						}
-                        contr = sld_in[i]-slope[i]*r;
-                        for (j=0; j<2; j++){
-                            if ( i == 0 && j == 0){
-                                continue;
-                                }
-                            if (i == n+1 && j == 1){
-                                continue;
-                                }
-                            if ( j == 1){
-                            	sign = 1.0;
-                                r += thick[i];
-                                }
-                            else{
-                            	sign = -1.0;
-								}
-
-                            qr = q * r;
-                            fun = 0.0;
-                            if(qr == 0.0){
-                                bes = sign * 1.0;
-                                }
-                            else{
-                                bes = sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
-                                if (fabs(slope[i]) > 0.0 ){
-                                    fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr);
-                                    }
-                                }
-                            vol = 4.0 * pi / 3.0 * r * r * r;
-                            //if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
-                            //	vol_sub += (vol_pre - vol);
-                            //}
-                            f += vol * (bes * contr + fun * slope[i]);
-                            }
-                            break;
-				}
-
-        }
-    //vol += vol_sub;
-    f2 = f * f / vol * 1.0e8;
-	f2 *= scale;
-	f2 += background;
-
-  free(sld_out);
-  free(slope);
-  free(sld_in);
-  free(thick);
-  free(A);
-  free(fun_type);
-
-    return (f2);
-}
-/**
- * Function to evaluate 1D scattering function
- * @param pars: parameters of the sphere
- * @param q: q-value
- * @return: function value
- */
-double onion_analytical_1D(OnionParameters *pars, double q) {
-	double dp[56];
-
-	dp[0] = pars->n_shells;
-	dp[1] = pars->scale;
-	dp[2] = pars->rad_core0;
-	dp[3] = pars->sld_core0;
-	dp[4] = pars->sld_solv;
-	dp[5] = pars->background;
-
-	dp[6] = pars->sld_out_shell1;
-	dp[7] = pars->sld_out_shell2;
-	dp[8] = pars->sld_out_shell3;
-	dp[9] = pars->sld_out_shell4;
-	dp[10] = pars->sld_out_shell5;
-	dp[11] = pars->sld_out_shell6;
-	dp[12] = pars->sld_out_shell7;
-	dp[13] = pars->sld_out_shell8;
-	dp[14] = pars->sld_out_shell9;
-	dp[15] = pars->sld_out_shell10;
-
-	dp[16] = pars->sld_in_shell1;
-	dp[17] = pars->sld_in_shell2;
-	dp[18] = pars->sld_in_shell3;
-	dp[19] = pars->sld_in_shell4;
-	dp[20] = pars->sld_in_shell5;
-	dp[21] = pars->sld_in_shell6;
-	dp[22] = pars->sld_in_shell7;
-	dp[23] = pars->sld_in_shell8;
-	dp[24] = pars->sld_in_shell9;
-	dp[25] = pars->sld_in_shell10;
-
-	dp[26] = pars->A_shell1;
-	dp[27] = pars->A_shell2;
-	dp[28] = pars->A_shell3;
-	dp[29] = pars->A_shell4;
-	dp[30] = pars->A_shell5;
-	dp[31] = pars->A_shell6;
-	dp[32] = pars->A_shell7;
-	dp[33] = pars->A_shell8;
-	dp[34] = pars->A_shell9;
-	dp[35] = pars->A_shell10;
-
-	dp[36] = pars->thick_shell1;
-	dp[37] = pars->thick_shell2;
-	dp[38] = pars->thick_shell3;
-	dp[39] = pars->thick_shell4;
-	dp[40] = pars->thick_shell5;
-	dp[41] = pars->thick_shell6;
-	dp[42] = pars->thick_shell7;
-	dp[43] = pars->thick_shell8;
-	dp[44] = pars->thick_shell9;
-	dp[45] = pars->thick_shell10;
-
-	dp[46] = pars->func_shell1;
-	dp[47] = pars->func_shell2;
-	dp[48] = pars->func_shell3;
-	dp[49] = pars->func_shell4;
-	dp[50] = pars->func_shell5;
-	dp[51] = pars->func_shell6;
-	dp[52] = pars->func_shell7;
-	dp[53] = pars->func_shell8;
-	dp[54] = pars->func_shell9;
-	dp[55] = pars->func_shell10;
-
-	return so_kernel(dp, q);
-}
-
-/**
- * Function to evaluate 2D scattering function
- * @param pars: parameters of the sphere
- * @param q: q-value
- * @return: function value
- */
-double onion_analytical_2D(OnionParameters *pars, double q, double phi) {
-	return onion_analytical_1D(pars,q);
-}
-
-double onion_analytical_2DXY(OnionParameters *pars, double qx, double qy){
-	return onion_analytical_1D(pars,sqrt(qx*qx+qy*qy));
-}
Index: sansmodels/src/c_extensions/onion.h
===================================================================
--- sansmodels/src/c_extensions/onion.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/onion.h	(revision c637521b288a6e3976fc9c0d461b44f70bb58c10)
@@ -1,4 +1,5 @@
 #if !defined(o_h)
 #define onion_h
+#include "parameters.hh"
 
 /**
@@ -34,138 +35,140 @@
  //[ORIENTATION_PARAMS]= <text> </text>
 
-typedef struct {
-	/// number of shells
-	//  [DEFAULT]=n_shells=1
-	int n_shells;
-    /// Scale factor
-    //  [DEFAULT]=scale= 1.0
-	double scale;
-    ///	Radius of sphere [A]
-    //  [DEFAULT]=rad_core0=200.0 [A]
-	double rad_core0;
-	///	sld_core [1/A^(2)]
-	//  [DEFAULT]=sld_core0= 1.0e-6 [1/A^(2)]
-	double sld_core0;
-	///	sld_solv [1/A^(2)]
-	//  [DEFAULT]=sld_solv= 6.4e-6 [1/A^(2)]
-	double sld_solv;
-	/// Incoherent Background [1/cm]
-	//  [DEFAULT]=background=0 [1/cm]
-	double background;
+class OnionModel{
+public:
+    // Model parameters
+    /// number of shells
+    //  [DEFAULT]=n_shells=1
+    Parameter n_shells;
+      /// Scale factor
+      //  [DEFAULT]=scale= 1.0
+    Parameter scale;
+      /// Radius of sphere [A]
+      //  [DEFAULT]=rad_core0=200.0 [A]
+    Parameter rad_core0;
+    /// sld_core [1/A^(2)]
+    //  [DEFAULT]=sld_core0= 1.0e-6 [1/A^(2)]
+    Parameter sld_core0;
+    /// sld_solv [1/A^(2)]
+    //  [DEFAULT]=sld_solv= 6.4e-6 [1/A^(2)]
+    Parameter sld_solv;
+    /// Incoherent Background [1/cm]
+    //  [DEFAULT]=background=0 [1/cm]
+    Parameter background;
 
     //  [DEFAULT]=sld_out_shell1=2.0e-06 [1/A^(2)]
-    double sld_out_shell1;
+    Parameter sld_out_shell1;
     //  [DEFAULT]=sld_out_shell2=2.5e-06 [1/A^(2)]
-    double sld_out_shell2;
+    Parameter sld_out_shell2;
     //  [DEFAULT]=sld_out_shell3=3.0e-06 [1/A^(2)]
-    double sld_out_shell3;
+    Parameter sld_out_shell3;
     //  [DEFAULT]=sld_out_shell4=3.5e-06 [1/A^(2)]
-    double sld_out_shell4;
+    Parameter sld_out_shell4;
     //  [DEFAULT]=sld_out_shell5=4.0e-06 [1/A^(2)]
-    double sld_out_shell5;
+    Parameter sld_out_shell5;
     //  [DEFAULT]=sld_out_shell6=4.5e-06 [1/A^(2)]
-    double sld_out_shell6;
+    Parameter sld_out_shell6;
     //  [DEFAULT]=sld_out_shell7=5.0e-06 [1/A^(2)]
-    double sld_out_shell7;
+    Parameter sld_out_shell7;
     //  [DEFAULT]=sld_out_shell8=5.5e-06 [1/A^(2)]
-    double sld_out_shell8;
+    Parameter sld_out_shell8;
     //  [DEFAULT]=sld_out_shell9=6.0e-06 [1/A^(2)]
-    double sld_out_shell9;
+    Parameter sld_out_shell9;
     //  [DEFAULT]=sld_out_shell10=6.2e-06 [1/A^(2)]
-    double sld_out_shell10;
+    Parameter sld_out_shell10;
 
     //  [DEFAULT]=sld_in_shell1=1.7e-06 [1/A^(2)]
-    double sld_in_shell1;
+    Parameter sld_in_shell1;
     //  [DEFAULT]=sld_in_shell2=2.2e-06 [1/A^(2)]
-    double sld_in_shell2;
+    Parameter sld_in_shell2;
     //  [DEFAULT]=sld_in_shell3=2.7e-06 [1/A^(2)]
-    double sld_in_shell3;
+    Parameter sld_in_shell3;
     //  [DEFAULT]=sld_in_shell4=3.2e-06 [1/A^(2)]
-    double sld_in_shell4;
+    Parameter sld_in_shell4;
     //  [DEFAULT]=sld_in_shell5=3.7e-06 [1/A^(2)]
-    double sld_in_shell5;
+    Parameter sld_in_shell5;
     //  [DEFAULT]=sld_in_shell6=4.2e-06 [1/A^(2)]
-    double sld_in_shell6;
+    Parameter sld_in_shell6;
     //  [DEFAULT]=sld_in_shell7=4.7e-06 [1/A^(2)]
-    double sld_in_shell7;
+    Parameter sld_in_shell7;
     //  [DEFAULT]=sld_in_shell8=5.2e-06 [1/A^(2)]
-    double sld_in_shell8;
+    Parameter sld_in_shell8;
     //  [DEFAULT]=sld_in_shell9=5.7e-06 [1/A^(2)]
-    double sld_in_shell9;
+    Parameter sld_in_shell9;
     //  [DEFAULT]=sld_in_shell10=6.0e-06 [1/A^(2)]
-    double sld_in_shell10;
+    Parameter sld_in_shell10;
 
     //  [DEFAULT]=A_shell1=1.0
-    double A_shell1;
+    Parameter A_shell1;
     //  [DEFAULT]=A_shell2=1.0
-    double A_shell2;
+    Parameter A_shell2;
     //  [DEFAULT]=A_shell3=1.0
-    double A_shell3;
+    Parameter A_shell3;
     //  [DEFAULT]=A_shell4=1.0
-    double A_shell4;
+    Parameter A_shell4;
     //  [DEFAULT]=A_shell5=1.0
-    double A_shell5;
+    Parameter A_shell5;
     //  [DEFAULT]=A_shell6=1.0
-    double A_shell6;
+    Parameter A_shell6;
     //  [DEFAULT]=A_shell7=1.0
-    double A_shell7;
+    Parameter A_shell7;
     //  [DEFAULT]=A_shell8=1.0
-    double A_shell8;
+    Parameter A_shell8;
     //  [DEFAULT]=A_shell9=1.0
-    double A_shell9;
+    Parameter A_shell9;
     //  [DEFAULT]=A_shell10=1.0
-    double A_shell10;
+    Parameter A_shell10;
 
     //  [DEFAULT]=thick_shell1=50.0 [A]
-    double thick_shell1;
+    Parameter thick_shell1;
     //  [DEFAULT]=thick_shell2=50.0 [A]
-    double thick_shell2;
+    Parameter thick_shell2;
     //  [DEFAULT]=thick_shell3=50.0 [A]
-    double thick_shell3;
+    Parameter thick_shell3;
     //  [DEFAULT]=thick_shell4=50.0 [A]
-    double thick_shell4;
+    Parameter thick_shell4;
     //  [DEFAULT]=thick_shell5=50.0 [A]
-    double thick_shell5;
+    Parameter thick_shell5;
     //  [DEFAULT]=thick_shell6=50.0 [A]
-    double thick_shell6;
+    Parameter thick_shell6;
     //  [DEFAULT]=thick_shell7=50.0 [A]
-    double thick_shell7;
+    Parameter thick_shell7;
     //  [DEFAULT]=thick_shell8=50.0 [A]
-    double thick_shell8;
+    Parameter thick_shell8;
     //  [DEFAULT]=thick_shell9=50.0 [A]
-    double thick_shell9;
+    Parameter thick_shell9;
     //  [DEFAULT]=thick_shell10=50.0 [A]
-    double thick_shell10;
+    Parameter thick_shell10;
 
     //  [DEFAULT]=func_shell1=2
-    int func_shell1;
+    Parameter func_shell1;
     //  [DEFAULT]=func_shell2=2
-    int func_shell2;
+    Parameter func_shell2;
     //  [DEFAULT]=func_shell3=2
-    int func_shell3;
+    Parameter func_shell3;
     //  [DEFAULT]=func_shell4=2
-    int func_shell4;
+    Parameter func_shell4;
     //  [DEFAULT]=func_shell5=2
-    int func_shell5;
+    Parameter func_shell5;
     //  [DEFAULT]=func_shell6=2
-    int func_shell6;
+    Parameter func_shell6;
     //  [DEFAULT]=func_shell7=2
-    int func_shell7;
+    Parameter func_shell7;
     //  [DEFAULT]=func_shell8=2
-    int func_shell8;
+    Parameter func_shell8;
     //  [DEFAULT]=func_shell9=2
-    int func_shell9;
+    Parameter func_shell9;
     //  [DEFAULT]=func_shell10=2
-    int func_shell10;
-} OnionParameters;
+    Parameter func_shell10;
 
-double so_kernel(double dq[], double q);
+  // Constructor
+  OnionModel();
 
-/// 1D scattering function
-double onion_analytical_1D(OnionParameters *pars, double q);
-
-/// 2D scattering function
-double onion_analytical_2D(OnionParameters *pars, double q, double phi);
-double onion_analytical_2DXY(OnionParameters *pars, double qx, double qy);
+  // 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/csparallelepiped.cpp
===================================================================
--- sansmodels/src/c_models/csparallelepiped.cpp	(revision df88829ada69c0e2705c41d796e0095d12a231db)
+++ sansmodels/src/c_models/csparallelepiped.cpp	(revision c637521b288a6e3976fc9c0d461b44f70bb58c10)
@@ -20,5 +20,4 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
@@ -26,32 +25,32 @@
 
 extern "C" {
-	#include "libCylinder.h"
-	#include "libStructureFactor.h"
-	#include "csparallelepiped.h"
-}
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+#include "csparallelepiped.h"
 
 // Convenience parameter structure
 typedef struct {
-    double scale;
-    double shortA;
-    double midB;
-    double longC;
-    double rimA;
-    double rimB;
-    double rimC;
-    double sld_rimA;
-    double sld_rimB;
-    double sld_rimC;
-    double sld_pcore;
-    double sld_solv;
-    double background;
-    double parallel_theta;
-    double parallel_phi;
-    double parallel_psi;
+  double scale;
+  double shortA;
+  double midB;
+  double longC;
+  double rimA;
+  double rimB;
+  double rimC;
+  double sld_rimA;
+  double sld_rimB;
+  double sld_rimC;
+  double sld_pcore;
+  double sld_solv;
+  double background;
+  double parallel_theta;
+  double parallel_phi;
+  double parallel_psi;
 } CSParallelepipedParameters;
 
 static double cspkernel(double dp[],double q, double ala, double alb, double alc){
-    // mu passed in is really mu*sqrt(1-sig^2)
-    double argA,argB,argC,argtA,argtB,argtC,tmp1,tmp2,tmp3,tmpt1,tmpt2,tmpt3;     //local variables
+  // mu passed in is really mu*sqrt(1-sig^2)
+  double argA,argB,argC,argtA,argtB,argtC,tmp1,tmp2,tmp3,tmpt1,tmpt2,tmpt3;     //local variables
 
   double aa,bb,cc, ta,tb,tc;
@@ -85,50 +84,50 @@
   tb=(aa+2.0*tb);///bb;
   tc=(aa+2.0*tc);
-    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    argA = q*aa*ala/2.0;
-    argB = q*bb*alb/2.0;
-    argC = q*cc*alc/2.0;
-    argtA = q*ta*ala/2.0;
+  //handle arg=0 separately, as sin(t)/t -> 1 as t->0
+  argA = q*aa*ala/2.0;
+  argB = q*bb*alb/2.0;
+  argC = q*cc*alc/2.0;
+  argtA = q*ta*ala/2.0;
   argtB = q*tb*alb/2.0;
   argtC = q*tc*alc/2.0;
 
-    if(argA==0.0) {
+  if(argA==0.0) {
     tmp1 = 1.0;
   } else {
     tmp1 = sin(argA)/argA;
-    }
-    if (argB==0.0) {
+  }
+  if (argB==0.0) {
     tmp2 = 1.0;
   } else {
     tmp2 = sin(argB)/argB;
-    }
-
-    if (argC==0.0) {
+  }
+
+  if (argC==0.0) {
     tmp3 = 1.0;
   } else {
     tmp3 = sin(argC)/argC;
-    }
-    if(argtA==0.0) {
+  }
+  if(argtA==0.0) {
     tmpt1 = 1.0;
   } else {
     tmpt1 = sin(argtA)/argtA;
-    }
-    if (argtB==0.0) {
+  }
+  if (argtB==0.0) {
     tmpt2 = 1.0;
   } else {
     tmpt2 = sin(argtB)/argtB;
-    }
-    if (argtC==0.0) {
+  }
+  if (argtC==0.0) {
     tmpt3 = 1.0;
   } else {
     tmpt3 = sin(argtC)*sin(argtC)/argtC/argtC;
-    }
-    // This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
-    retVal =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)*
-        ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
-    //retVal *= (tmp3*tmp3);
-    retVal /= Vot;
-
-    return (retVal);
+  }
+  // This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
+  retVal =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)*
+      ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
+  //retVal *= (tmp3*tmp3);
+  retVal /= Vot;
+
+  return (retVal);
 
 }//Function cspkernel()
@@ -176,38 +175,38 @@
 
 
-    // parallelepiped c axis orientation
-    cparallel_x = sin(theta) * cos(phi);
-    cparallel_y = sin(theta) * sin(phi);
-    cparallel_z = cos(theta);
-
-    // q vector
-    q_z = 0.0;
-
-    // Compute the angle btw vector q and the
-    // axis of the parallelepiped
-    cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z;
-    alpha = acos(cos_val_c);
-
-    // parallelepiped a axis orientation
-    parallel_x = sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi);
-    parallel_y = cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi);
-
-    cos_val_a = parallel_x*q_x + parallel_y*q_y;
-
-
-
-    // parallelepiped b axis orientation
-    bparallel_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi);
-    bparallel_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi);
-    // axis of the parallelepiped
-    cos_val_b = sin(acos(cos_val_a)) ;
-
-
-
-    // The following test should always pass
-    if (fabs(cos_val_c)>1.0) {
-      printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
-      return 0;
-    }
+  // parallelepiped c axis orientation
+  cparallel_x = sin(theta) * cos(phi);
+  cparallel_y = sin(theta) * sin(phi);
+  cparallel_z = cos(theta);
+
+  // q vector
+  q_z = 0.0;
+
+  // Compute the angle btw vector q and the
+  // axis of the parallelepiped
+  cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z;
+  alpha = acos(cos_val_c);
+
+  // parallelepiped a axis orientation
+  parallel_x = sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi);
+  parallel_y = cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi);
+
+  cos_val_a = parallel_x*q_x + parallel_y*q_y;
+
+
+
+  // parallelepiped b axis orientation
+  bparallel_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi);
+  bparallel_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi);
+  // axis of the parallelepiped
+  cos_val_b = sin(acos(cos_val_a)) ;
+
+
+
+  // The following test should always pass
+  if (fabs(cos_val_c)>1.0) {
+    printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
 
   // Call the IGOR library function to get the kernel
@@ -235,5 +234,5 @@
   double q;
   q = sqrt(qx*qx+qy*qy);
-    return csparallelepiped_analytical_2D_scaled(pars, q, qx/q, qy/q);
+  return csparallelepiped_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
@@ -242,23 +241,23 @@
 
 CSParallelepipedModel :: CSParallelepipedModel() {
-	scale      = Parameter(1.0);
-	shortA     = Parameter(35.0, true);
-	shortA.set_min(1.0);
-	midB     = Parameter(75.0, true);
-	midB.set_min(1.0);
-	longC    = Parameter(400.0, true);
-	longC.set_min(1.0);
-	rimA     = Parameter(10.0, true);
-	rimB     = Parameter(10.0, true);
-	rimC     = Parameter(10.0, true);
-	sld_rimA     = Parameter(2.0e-6, true);
-	sld_rimB     = Parameter(4.0e-6, true);
-	sld_rimC    = Parameter(2.0e-6, true);
-	sld_pcore   = Parameter(1.0e-6);
-	sld_solv   = Parameter(6.0e-6);
-	background = Parameter(0.06);
-	parallel_theta  = Parameter(0.0, true);
-	parallel_phi    = Parameter(0.0, true);
-	parallel_psi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  shortA     = Parameter(35.0, true);
+  shortA.set_min(1.0);
+  midB     = Parameter(75.0, true);
+  midB.set_min(1.0);
+  longC    = Parameter(400.0, true);
+  longC.set_min(1.0);
+  rimA     = Parameter(10.0, true);
+  rimB     = Parameter(10.0, true);
+  rimC     = Parameter(10.0, true);
+  sld_rimA     = Parameter(2.0e-6, true);
+  sld_rimB     = Parameter(4.0e-6, true);
+  sld_rimC    = Parameter(2.0e-6, true);
+  sld_pcore   = Parameter(1.0e-6);
+  sld_solv   = Parameter(6.0e-6);
+  background = Parameter(0.06);
+  parallel_theta  = Parameter(0.0, true);
+  parallel_phi    = Parameter(0.0, true);
+  parallel_psi    = Parameter(0.0, true);
 }
 
@@ -270,73 +269,73 @@
  */
 double CSParallelepipedModel :: operator()(double q) {
-	double dp[13];
-
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = shortA();
-	dp[2] = midB();
-	dp[3] = longC();
-	dp[4] = rimA();
-	dp[5] = rimB();
-	dp[6] = rimC();
-	dp[7] = sld_rimA();
-	dp[8] = sld_rimB();
-	dp[9] = sld_rimC();
-	dp[10] = sld_pcore();
-	dp[11] = sld_solv();
-	dp[12] = 0.0;
-
-	// Get the dispersion points for the short_edgeA
-	vector<WeightPoint> weights_shortA;
-	shortA.get_weights(weights_shortA);
-
-	// Get the dispersion points for the longer_edgeB
-	vector<WeightPoint> weights_midB;
-	midB.get_weights(weights_midB);
-
-	// Get the dispersion points for the longuest_edgeC
-	vector<WeightPoint> weights_longC;
-	longC.get_weights(weights_longC);
-
-
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over short_edgeA weight points
-	for(int i=0; i< (int)weights_shortA.size(); i++) {
-		dp[1] = weights_shortA[i].value;
-
-		// Loop over longer_edgeB weight points
-		for(int j=0; j< (int)weights_midB.size(); j++) {
-			dp[2] = weights_midB[j].value;
-
-			// Loop over longuest_edgeC weight points
-			for(int k=0; k< (int)weights_longC.size(); k++) {
-				dp[3] = weights_longC[k].value;
-				//Un-normalize  by volume
-				sum += weights_shortA[i].weight * weights_midB[j].weight
-					* weights_longC[k].weight * CSParallelepiped(dp, q)
-					* weights_shortA[i].value*weights_midB[j].value
-					* weights_longC[k].value;
-				//Find average volume
-				vol += weights_shortA[i].weight * weights_midB[j].weight
-					* weights_longC[k].weight
-					* weights_shortA[i].value * weights_midB[j].value
-					* weights_longC[k].value;
-
-				norm += weights_shortA[i].weight
-					 * weights_midB[j].weight * weights_longC[k].weight;
-			}
-		}
-	}
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[13];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = shortA();
+  dp[2] = midB();
+  dp[3] = longC();
+  dp[4] = rimA();
+  dp[5] = rimB();
+  dp[6] = rimC();
+  dp[7] = sld_rimA();
+  dp[8] = sld_rimB();
+  dp[9] = sld_rimC();
+  dp[10] = sld_pcore();
+  dp[11] = sld_solv();
+  dp[12] = 0.0;
+
+  // Get the dispersion points for the short_edgeA
+  vector<WeightPoint> weights_shortA;
+  shortA.get_weights(weights_shortA);
+
+  // Get the dispersion points for the longer_edgeB
+  vector<WeightPoint> weights_midB;
+  midB.get_weights(weights_midB);
+
+  // Get the dispersion points for the longuest_edgeC
+  vector<WeightPoint> weights_longC;
+  longC.get_weights(weights_longC);
+
+
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over short_edgeA weight points
+  for(int i=0; i< (int)weights_shortA.size(); i++) {
+    dp[1] = weights_shortA[i].value;
+
+    // Loop over longer_edgeB weight points
+    for(int j=0; j< (int)weights_midB.size(); j++) {
+      dp[2] = weights_midB[j].value;
+
+      // Loop over longuest_edgeC weight points
+      for(int k=0; k< (int)weights_longC.size(); k++) {
+        dp[3] = weights_longC[k].value;
+        //Un-normalize  by volume
+        sum += weights_shortA[i].weight * weights_midB[j].weight
+            * weights_longC[k].weight * CSParallelepiped(dp, q)
+        * weights_shortA[i].value*weights_midB[j].value
+        * weights_longC[k].value;
+        //Find average volume
+        vol += weights_shortA[i].weight * weights_midB[j].weight
+            * weights_longC[k].weight
+            * weights_shortA[i].value * weights_midB[j].value
+            * weights_longC[k].value;
+
+        norm += weights_shortA[i].weight
+            * weights_midB[j].weight * weights_longC[k].weight;
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 /**
@@ -347,129 +346,129 @@
  */
 double CSParallelepipedModel :: operator()(double qx, double qy) {
-	CSParallelepipedParameters dp;
-	// Fill parameter array
-	dp.scale      = scale();
-	dp.shortA   = shortA();
-	dp.midB   = midB();
-	dp.longC  = longC();
-	dp.rimA   = rimA();
-	dp.rimB   = rimB();
-	dp.rimC  = rimC();
-	dp.sld_rimA   = sld_rimA();
-	dp.sld_rimB   = sld_rimB();
-	dp.sld_rimC  = sld_rimC();
-	dp.sld_pcore   = sld_pcore();
-	dp.sld_solv   = sld_solv();
-	dp.background = 0.0;
-	//dp.background = background();
-	dp.parallel_theta  = parallel_theta();
-	dp.parallel_phi    = parallel_phi();
-	dp.parallel_psi    = parallel_psi();
-
-
-
-	// Get the dispersion points for the short_edgeA
-	vector<WeightPoint> weights_shortA;
-	shortA.get_weights(weights_shortA);
-
-	// Get the dispersion points for the longer_edgeB
-	vector<WeightPoint> weights_midB;
-	midB.get_weights(weights_midB);
-
-	// Get the dispersion points for the longuest_edgeC
-	vector<WeightPoint> weights_longC;
-	longC.get_weights(weights_longC);
-
-	// Get angular averaging for theta
-	vector<WeightPoint> weights_parallel_theta;
-	parallel_theta.get_weights(weights_parallel_theta);
-
-	// Get angular averaging for phi
-	vector<WeightPoint> weights_parallel_phi;
-	parallel_phi.get_weights(weights_parallel_phi);
-
-	// Get angular averaging for psi
-	vector<WeightPoint> weights_parallel_psi;
-	parallel_psi.get_weights(weights_parallel_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(int i=0; i< (int)weights_shortA.size(); i++) {
-		dp.shortA = weights_shortA[i].value;
-
-		// Loop over longer_edgeB weight points
-		for(int j=0; j< (int)weights_midB.size(); j++) {
-			dp.midB = weights_midB[j].value;
-
-			// Average over longuest_edgeC distribution
-			for(int k=0; k< (int)weights_longC.size(); k++) {
-				dp.longC = weights_longC[k].value;
-
-				// Average over theta distribution
-				for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
-				dp.parallel_theta = weights_parallel_theta[l].value;
-
-					// Average over phi distribution
-					for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
-						dp.parallel_phi = weights_parallel_phi[m].value;
-
-						// Average over phi distribution
-						for(int n=0; n< (int)weights_parallel_psi.size(); n++) {
-							dp.parallel_psi = weights_parallel_psi[n].value;
-							//Un-normalize by volume
-							double _ptvalue = weights_shortA[i].weight
-								* weights_midB[j].weight
-								* weights_longC[k].weight
-								* weights_parallel_theta[l].weight
-								* weights_parallel_phi[m].weight
-								* weights_parallel_psi[n].weight
-								* csparallelepiped_analytical_2DXY(&dp, qx, qy)
-								* weights_shortA[i].value*weights_midB[j].value
-								* weights_longC[k].value;
-
-							if (weights_parallel_theta.size()>1) {
-								_ptvalue *= fabs(sin(weights_parallel_theta[l].value*pi/180.0));
-							}
-							sum += _ptvalue;
-							//Find average volume
-							vol += weights_shortA[i].weight
-								* weights_midB[j].weight
-								* weights_longC[k].weight
-								* weights_shortA[i].value*weights_midB[j].value
-								* weights_longC[k].value;
-							//Find norm for volume
-							norm_vol += weights_shortA[i].weight
-								* weights_midB[j].weight
-								* weights_longC[k].weight;
-
-							norm += weights_shortA[i].weight
-								* weights_midB[j].weight
-								* weights_longC[k].weight
-								* weights_parallel_theta[l].weight
-								* weights_parallel_phi[m].weight
-								* weights_parallel_psi[n].weight;
-						}
-					}
-
-				}
-			}
-		}
-	}
-	// Averaging in theta needs an extra normalization
-	// factor to account for the sin(theta) term in the
-	// integration (see documentation).
-	if (weights_parallel_theta.size()>1) norm = norm / asin(1.0);
-
-	if (vol != 0.0 && norm_vol != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm_vol);}
-
-	return sum/norm + background();
+  CSParallelepipedParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.shortA   = shortA();
+  dp.midB   = midB();
+  dp.longC  = longC();
+  dp.rimA   = rimA();
+  dp.rimB   = rimB();
+  dp.rimC  = rimC();
+  dp.sld_rimA   = sld_rimA();
+  dp.sld_rimB   = sld_rimB();
+  dp.sld_rimC  = sld_rimC();
+  dp.sld_pcore   = sld_pcore();
+  dp.sld_solv   = sld_solv();
+  dp.background = 0.0;
+  //dp.background = background();
+  dp.parallel_theta  = parallel_theta();
+  dp.parallel_phi    = parallel_phi();
+  dp.parallel_psi    = parallel_psi();
+
+
+
+  // Get the dispersion points for the short_edgeA
+  vector<WeightPoint> weights_shortA;
+  shortA.get_weights(weights_shortA);
+
+  // Get the dispersion points for the longer_edgeB
+  vector<WeightPoint> weights_midB;
+  midB.get_weights(weights_midB);
+
+  // Get the dispersion points for the longuest_edgeC
+  vector<WeightPoint> weights_longC;
+  longC.get_weights(weights_longC);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_parallel_theta;
+  parallel_theta.get_weights(weights_parallel_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_parallel_phi;
+  parallel_phi.get_weights(weights_parallel_phi);
+
+  // Get angular averaging for psi
+  vector<WeightPoint> weights_parallel_psi;
+  parallel_psi.get_weights(weights_parallel_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(int i=0; i< (int)weights_shortA.size(); i++) {
+    dp.shortA = weights_shortA[i].value;
+
+    // Loop over longer_edgeB weight points
+    for(int j=0; j< (int)weights_midB.size(); j++) {
+      dp.midB = weights_midB[j].value;
+
+      // Average over longuest_edgeC distribution
+      for(int k=0; k< (int)weights_longC.size(); k++) {
+        dp.longC = weights_longC[k].value;
+
+        // Average over theta distribution
+        for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
+          dp.parallel_theta = weights_parallel_theta[l].value;
+
+          // Average over phi distribution
+          for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
+            dp.parallel_phi = weights_parallel_phi[m].value;
+
+            // Average over phi distribution
+            for(int n=0; n< (int)weights_parallel_psi.size(); n++) {
+              dp.parallel_psi = weights_parallel_psi[n].value;
+              //Un-normalize by volume
+              double _ptvalue = weights_shortA[i].weight
+                  * weights_midB[j].weight
+                  * weights_longC[k].weight
+                  * weights_parallel_theta[l].weight
+                  * weights_parallel_phi[m].weight
+                  * weights_parallel_psi[n].weight
+                  * csparallelepiped_analytical_2DXY(&dp, qx, qy)
+              * weights_shortA[i].value*weights_midB[j].value
+              * weights_longC[k].value;
+
+              if (weights_parallel_theta.size()>1) {
+                _ptvalue *= fabs(sin(weights_parallel_theta[l].value*pi/180.0));
+              }
+              sum += _ptvalue;
+              //Find average volume
+              vol += weights_shortA[i].weight
+                  * weights_midB[j].weight
+                  * weights_longC[k].weight
+                  * weights_shortA[i].value*weights_midB[j].value
+                  * weights_longC[k].value;
+              //Find norm for volume
+              norm_vol += weights_shortA[i].weight
+                  * weights_midB[j].weight
+                  * weights_longC[k].weight;
+
+              norm += weights_shortA[i].weight
+                  * weights_midB[j].weight
+                  * weights_longC[k].weight
+                  * weights_parallel_theta[l].weight
+                  * weights_parallel_phi[m].weight
+                  * weights_parallel_psi[n].weight;
+            }
+          }
+
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_parallel_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -483,7 +482,7 @@
  */
 double CSParallelepipedModel :: evaluate_rphi(double q, double phi) {
-	double qx = q*cos(phi);
-	double qy = q*sin(phi);
-	return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -492,63 +491,63 @@
  */
 double CSParallelepipedModel :: calculate_ER() {
-	CSParallelepipedParameters dp;
-	dp.shortA   = shortA();
-	dp.midB   = midB();
-	dp.longC  = longC();
-	dp.rimA   = rimA();
-	dp.rimB   = rimB();
-	dp.rimC  = rimC();
-
-	double rad_out = 0.0;
-	double pi = 4.0*atan(1.0);
-	double suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
-	double height =(dp.longC + 2.0*dp.rimC);
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the short_edgeA
-	vector<WeightPoint> weights_shortA;
-	shortA.get_weights(weights_shortA);
-
-	// Get the dispersion points for the longer_edgeB
-	vector<WeightPoint> weights_midB;
-	midB.get_weights(weights_midB);
-
-	// Get angular averaging for the longuest_edgeC
-	vector<WeightPoint> weights_longC;
-	longC.get_weights(weights_longC);
-
-	// Loop over radius weight points
-	for(int i=0; i< (int)weights_shortA.size(); i++) {
-		dp.shortA = weights_shortA[i].value;
-
-		// Loop over longer_edgeB weight points
-		for(int j=0; j< (int)weights_midB.size(); j++) {
-			dp.midB = weights_midB[j].value;
-
-			// Average over longuest_edgeC distribution
-			for(int k=0; k< (int)weights_longC.size(); k++) {
-				dp.longC = weights_longC[k].value;
-				//Calculate surface averaged radius
-				//This is rough approximation.
-				suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
-				height =(dp.longC + 2.0*dp.rimC);
-				//Note: output of "DiamCyl(dp.length,dp.radius)" is a DIAMETER.
-				sum +=weights_shortA[i].weight* weights_midB[j].weight
-					* weights_longC[k].weight*DiamCyl(height, suf_rad)/2.0;
-				norm += weights_shortA[i].weight* weights_midB[j].weight*weights_longC[k].weight;
-			}
-		}
-	}
-
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		//Note: output of "DiamCyl(length,radius)" is DIAMETER.
-		rad_out = DiamCyl(dp.longC, suf_rad)/2.0;}
-	return rad_out;
-
-}
+  CSParallelepipedParameters dp;
+  dp.shortA   = shortA();
+  dp.midB   = midB();
+  dp.longC  = longC();
+  dp.rimA   = rimA();
+  dp.rimB   = rimB();
+  dp.rimC  = rimC();
+
+  double rad_out = 0.0;
+  double pi = 4.0*atan(1.0);
+  double suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
+  double height =(dp.longC + 2.0*dp.rimC);
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the short_edgeA
+  vector<WeightPoint> weights_shortA;
+  shortA.get_weights(weights_shortA);
+
+  // Get the dispersion points for the longer_edgeB
+  vector<WeightPoint> weights_midB;
+  midB.get_weights(weights_midB);
+
+  // Get angular averaging for the longuest_edgeC
+  vector<WeightPoint> weights_longC;
+  longC.get_weights(weights_longC);
+
+  // Loop over radius weight points
+  for(int i=0; i< (int)weights_shortA.size(); i++) {
+    dp.shortA = weights_shortA[i].value;
+
+    // Loop over longer_edgeB weight points
+    for(int j=0; j< (int)weights_midB.size(); j++) {
+      dp.midB = weights_midB[j].value;
+
+      // Average over longuest_edgeC distribution
+      for(int k=0; k< (int)weights_longC.size(); k++) {
+        dp.longC = weights_longC[k].value;
+        //Calculate surface averaged radius
+        //This is rough approximation.
+        suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
+        height =(dp.longC + 2.0*dp.rimC);
+        //Note: output of "DiamCyl(dp.length,dp.radius)" is a DIAMETER.
+        sum +=weights_shortA[i].weight* weights_midB[j].weight
+            * weights_longC[k].weight*DiamCyl(height, suf_rad)/2.0;
+        norm += weights_shortA[i].weight* weights_midB[j].weight*weights_longC[k].weight;
+      }
+    }
+  }
+
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl(length,radius)" is DIAMETER.
+    rad_out = DiamCyl(dp.longC, suf_rad)/2.0;}
+  return rad_out;
+
+}
Index: sansmodels/src/c_models/models.hh
===================================================================
--- sansmodels/src/c_models/models.hh	(revision df88829ada69c0e2705c41d796e0095d12a231db)
+++ sansmodels/src/c_models/models.hh	(revision c637521b288a6e3976fc9c0d461b44f70bb58c10)
@@ -28,78 +28,5 @@
 
 
-class OnionModel{
-public:
-	// Model parameters
-	Parameter n_shells;
-	Parameter scale;
-	Parameter rad_core0;
-	Parameter sld_core0;
-	Parameter sld_solv;
-	Parameter background;
-
-	Parameter sld_out_shell1;
-	Parameter sld_out_shell2;
-	Parameter sld_out_shell3;
-	Parameter sld_out_shell4;
-	Parameter sld_out_shell5;
-	Parameter sld_out_shell6;
-	Parameter sld_out_shell7;
-	Parameter sld_out_shell8;
-	Parameter sld_out_shell9;
-	Parameter sld_out_shell10;
-
-	Parameter sld_in_shell1;
-	Parameter sld_in_shell2;
-	Parameter sld_in_shell3;
-	Parameter sld_in_shell4;
-	Parameter sld_in_shell5;
-	Parameter sld_in_shell6;
-	Parameter sld_in_shell7;
-	Parameter sld_in_shell8;
-	Parameter sld_in_shell9;
-	Parameter sld_in_shell10;
-
-	Parameter A_shell1;
-	Parameter A_shell2;
-	Parameter A_shell3;
-	Parameter A_shell4;
-	Parameter A_shell5;
-	Parameter A_shell6;
-	Parameter A_shell7;
-	Parameter A_shell8;
-	Parameter A_shell9;
-	Parameter A_shell10;
-
-	Parameter thick_shell1;
-	Parameter thick_shell2;
-	Parameter thick_shell3;
-	Parameter thick_shell4;
-	Parameter thick_shell5;
-	Parameter thick_shell6;
-	Parameter thick_shell7;
-	Parameter thick_shell8;
-	Parameter thick_shell9;
-	Parameter thick_shell10;
-
-	Parameter func_shell1;
-	Parameter func_shell2;
-	Parameter func_shell3;
-	Parameter func_shell4;
-	Parameter func_shell5;
-	Parameter func_shell6;
-	Parameter func_shell7;
-	Parameter func_shell8;
-	Parameter func_shell9;
-	Parameter func_shell10;
-
-	// Constructor
-	OnionModel();
-
-	// 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/onion.cpp
===================================================================
--- sansmodels/src/c_models/onion.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/onion.cpp	(revision c637521b288a6e3976fc9c0d461b44f70bb58c10)
@@ -21,81 +21,309 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
-
-extern "C" {
-	#include "onion.h"
+#include "onion.h"
+
+// Convenience parameter structure
+typedef struct {
+  int n_shells;
+  double scale;
+  double rad_core0;
+  double sld_core0;
+  double sld_solv;
+  double background;
+  double sld_out_shell1;
+  double sld_out_shell2;
+  double sld_out_shell3;
+  double sld_out_shell4;
+  double sld_out_shell5;
+  double sld_out_shell6;
+  double sld_out_shell7;
+  double sld_out_shell8;
+  double sld_out_shell9;
+  double sld_out_shell10;
+  double sld_in_shell1;
+  double sld_in_shell2;
+  double sld_in_shell3;
+  double sld_in_shell4;
+  double sld_in_shell5;
+  double sld_in_shell6;
+  double sld_in_shell7;
+  double sld_in_shell8;
+  double sld_in_shell9;
+  double sld_in_shell10;
+
+  double A_shell1;
+  double A_shell2;
+  double A_shell3;
+  double A_shell4;
+  double A_shell5;
+  double A_shell6;
+  double A_shell7;
+  double A_shell8;
+  double A_shell9;
+  double A_shell10;
+
+  double thick_shell1;
+  double thick_shell2;
+  double thick_shell3;
+  double thick_shell4;
+  double thick_shell5;
+  double thick_shell6;
+  double thick_shell7;
+  double thick_shell8;
+  double thick_shell9;
+  double thick_shell10;
+
+  int func_shell1;
+  int func_shell2;
+  int func_shell3;
+  int func_shell4;
+  int func_shell5;
+  int func_shell6;
+  int func_shell7;
+  int func_shell8;
+  int func_shell9;
+  int func_shell10;
+} OnionParameters;
+
+// some details can be found in sld_cal.c
+static double so_kernel(double dp[], double q) {
+  int n = dp[0];
+  double scale = dp[1];
+  double rad_core0 = dp[2];
+  double sld_core0 = dp[3];
+  double sld_solv = dp[4];
+  double background = dp[5];
+  int i,j;
+  double bes,fun,alpha,f,vol,vol_pre,vol_sub,qr,r,contr,f2;
+  double sign;
+  double pi;
+  double r0 = 0.0;
+
+  double *sld_out;
+  double *slope;
+  double *sld_in;
+  double *thick;
+  double *A;
+  int *fun_type;
+
+  sld_out = (double*)malloc((n+2)*sizeof(double));
+  slope = (double*)malloc((n+2)*sizeof(double));
+  sld_in = (double*)malloc((n+2)*sizeof(double));
+  thick = (double*)malloc((n+2)*sizeof(double));
+  A = (double*)malloc((n+2)*sizeof(double));
+  fun_type = (int*)malloc((n+2)*sizeof(int));
+
+  for (i =1; i<=n; i++){
+    sld_out[i] = dp[i+5];
+    sld_in[i] = dp[i+15];
+    A[i] = dp[i+25];
+    thick[i] = dp[i+35];
+    fun_type[i] = dp[i+45];
+  }
+  sld_out[0] = sld_core0;
+  sld_out[n+1] = sld_solv;
+  sld_in[0] = sld_core0;
+  sld_in[n+1] = sld_solv;
+  thick[0] = rad_core0;
+  thick[n+1] = 1e+10;
+  A[0] = 0.0;
+  A[n+1] = 0.0;
+  fun_type[0] = 0;
+  fun_type[n+1] = 0;
+
+
+  pi = 4.0*atan(1.0);
+  f = 0.0;
+  r = 0.0;
+  vol = 0.0;
+  vol_pre = 0.0;
+  vol_sub = 0.0;
+
+  for (i =0; i<= n+1; i++){
+    if (thick[i] == 0.0){
+      continue;
+    }
+    if (fun_type[i]== 0 ){
+      slope[i] = 0.0;
+      A[i] = 0.0;
+    }
+    vol_pre = vol;
+    switch(fun_type[i]){
+    case 2 :
+      r0 = r;
+      if (A[i] == 0.0){
+        slope[i] = 0.0;
+      }
+      else{
+        slope[i] = (sld_out[i]-sld_in[i])/(exp(A[i])-1.0);
+      }
+      for (j=0; j<2; j++){
+        if ( i == 0 && j == 0){
+          continue;
+        }
+        if (i == n+1 && j == 1){
+          continue;
+        }
+        if ( j == 1){
+          sign = 1.0;
+          r += thick[i];
+        }
+        else{
+          sign = -1.0;
+        }
+        qr = q * r;
+        alpha = A[i] * r/thick[i];
+        fun = 0.0;
+        if(qr == 0.0){
+          fun = sign * 1.0;
+          bes = sign * 1.0;
+        }
+        else{
+          if (fabs(A[i]) > 0.0 ){
+            fun = 3.0 * ((alpha*alpha - qr * qr) * sin(qr) - 2.0 * alpha * qr * cos(qr))/ ((alpha * alpha + qr * qr) * (alpha * alpha + qr * qr) * qr);
+            fun = fun - 3.0 * (alpha * sin(qr) - qr * cos(qr)) / ((alpha * alpha + qr * qr) * qr);
+            fun = - sign *fun;
+            bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
+          }
+          else {
+            fun = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
+            bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
+          }
+        }
+        contr = slope[i]*exp(A[i]*(r-r0)/thick[i]);
+
+        vol = 4.0 * pi / 3.0 * r * r * r;
+        //if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && A[i]==0.0){
+        //  vol_sub += (vol_pre - vol);
+        //}
+        f += vol * (contr * (fun) + (sld_in[i]-slope[i]) * bes);
+      }
+      break;
+    default :
+      if (fun_type[i]==0){
+        slope[i] = 0.0;
+      }
+      else{
+        slope[i]= (sld_out[i] -sld_in[i])/thick[i];
+      }
+      contr = sld_in[i]-slope[i]*r;
+      for (j=0; j<2; j++){
+        if ( i == 0 && j == 0){
+          continue;
+        }
+        if (i == n+1 && j == 1){
+          continue;
+        }
+        if ( j == 1){
+          sign = 1.0;
+          r += thick[i];
+        }
+        else{
+          sign = -1.0;
+        }
+
+        qr = q * r;
+        fun = 0.0;
+        if(qr == 0.0){
+          bes = sign * 1.0;
+        }
+        else{
+          bes = sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
+          if (fabs(slope[i]) > 0.0 ){
+            fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr);
+          }
+        }
+        vol = 4.0 * pi / 3.0 * r * r * r;
+        //if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
+        //  vol_sub += (vol_pre - vol);
+        //}
+        f += vol * (bes * contr + fun * slope[i]);
+      }
+      break;
+    }
+
+  }
+  //vol += vol_sub;
+  f2 = f * f / vol * 1.0e8;
+  f2 *= scale;
+  f2 += background;
+
+  free(sld_out);
+  free(slope);
+  free(sld_in);
+  free(thick);
+  free(A);
+  free(fun_type);
+
+  return (f2);
 }
 
+
 OnionModel :: OnionModel() {
-	n_shells = Parameter(1.0);
-	scale = Parameter(1.0);
-	rad_core0 = Parameter(200.0);
-	sld_core0 = Parameter(1e-06);
-	sld_solv = Parameter(6.4e-06);
-    background = Parameter(0.0);
-
-
-    sld_out_shell1 = Parameter(1.0e-06);
-    sld_out_shell2 = Parameter(1.0e-06);
-    sld_out_shell3 = Parameter(1.0e-06);
-    sld_out_shell4 = Parameter(1.0e-06);
-    sld_out_shell5 = Parameter(1.0e-06);
-    sld_out_shell6 = Parameter(1.0e-06);
-    sld_out_shell7 = Parameter(1.0e-06);
-    sld_out_shell8 = Parameter(1.0e-06);
-    sld_out_shell9 = Parameter(1.0e-06);
-    sld_out_shell10 = Parameter(1.0e-06);
-
-
-    sld_in_shell1 = Parameter(2.3e-06);
-    sld_in_shell2 = Parameter(2.6e-06);
-    sld_in_shell3 = Parameter(2.9e-06);
-    sld_in_shell4 = Parameter(3.2e-06);
-    sld_in_shell5 = Parameter(3.5e-06);
-    sld_in_shell6 = Parameter(3.8e-06);
-    sld_in_shell7 = Parameter(4.1e-06);
-    sld_in_shell8 = Parameter(4.4e-06);
-    sld_in_shell9 = Parameter(4.7e-06);
-    sld_in_shell10 = Parameter(5.0e-06);
-
-
-    A_shell1 = Parameter(1.0);
-    A_shell2 = Parameter(1.0);
-    A_shell3 = Parameter(1.0);
-    A_shell4 = Parameter(1.0);
-    A_shell5 = Parameter(1.0);
-    A_shell6 = Parameter(1.0);
-    A_shell7 = Parameter(1.0);
-    A_shell8 = Parameter(1.0);
-    A_shell9 = Parameter(1.0);
-    A_shell10 = Parameter(1.0);
-
-
-    thick_shell1 = Parameter(50.0);
-    thick_shell2 = Parameter(50.0);
-    thick_shell3 = Parameter(50.0);
-    thick_shell4 = Parameter(50.0);
-    thick_shell5 = Parameter(50.0);
-    thick_shell6 = Parameter(50.0);
-    thick_shell7 = Parameter(50.0);
-    thick_shell8 = Parameter(50.0);
-    thick_shell9 = Parameter(50.0);
-    thick_shell10 = Parameter(50.0);
-
-
-    func_shell1 = Parameter(2);
-    func_shell2 = Parameter(2);
-    func_shell3 = Parameter(2);
-    func_shell4 = Parameter(2);
-    func_shell5 = Parameter(2);
-    func_shell6 = Parameter(2);
-    func_shell7 = Parameter(2);
-    func_shell8 = Parameter(2);
-    func_shell9 = Parameter(2);
-    func_shell10 = Parameter(2);
-
+  n_shells = Parameter(1.0);
+  scale = Parameter(1.0);
+  rad_core0 = Parameter(200.0);
+  sld_core0 = Parameter(1e-06);
+  sld_solv = Parameter(6.4e-06);
+  background = Parameter(0.0);
+
+  sld_out_shell1 = Parameter(1.0e-06);
+  sld_out_shell2 = Parameter(1.0e-06);
+  sld_out_shell3 = Parameter(1.0e-06);
+  sld_out_shell4 = Parameter(1.0e-06);
+  sld_out_shell5 = Parameter(1.0e-06);
+  sld_out_shell6 = Parameter(1.0e-06);
+  sld_out_shell7 = Parameter(1.0e-06);
+  sld_out_shell8 = Parameter(1.0e-06);
+  sld_out_shell9 = Parameter(1.0e-06);
+  sld_out_shell10 = Parameter(1.0e-06);
+
+  sld_in_shell1 = Parameter(2.3e-06);
+  sld_in_shell2 = Parameter(2.6e-06);
+  sld_in_shell3 = Parameter(2.9e-06);
+  sld_in_shell4 = Parameter(3.2e-06);
+  sld_in_shell5 = Parameter(3.5e-06);
+  sld_in_shell6 = Parameter(3.8e-06);
+  sld_in_shell7 = Parameter(4.1e-06);
+  sld_in_shell8 = Parameter(4.4e-06);
+  sld_in_shell9 = Parameter(4.7e-06);
+  sld_in_shell10 = Parameter(5.0e-06);
+
+  A_shell1 = Parameter(1.0);
+  A_shell2 = Parameter(1.0);
+  A_shell3 = Parameter(1.0);
+  A_shell4 = Parameter(1.0);
+  A_shell5 = Parameter(1.0);
+  A_shell6 = Parameter(1.0);
+  A_shell7 = Parameter(1.0);
+  A_shell8 = Parameter(1.0);
+  A_shell9 = Parameter(1.0);
+  A_shell10 = Parameter(1.0);
+
+  thick_shell1 = Parameter(50.0);
+  thick_shell2 = Parameter(50.0);
+  thick_shell3 = Parameter(50.0);
+  thick_shell4 = Parameter(50.0);
+  thick_shell5 = Parameter(50.0);
+  thick_shell6 = Parameter(50.0);
+  thick_shell7 = Parameter(50.0);
+  thick_shell8 = Parameter(50.0);
+  thick_shell9 = Parameter(50.0);
+  thick_shell10 = Parameter(50.0);
+
+  func_shell1 = Parameter(2);
+  func_shell2 = Parameter(2);
+  func_shell3 = Parameter(2);
+  func_shell4 = Parameter(2);
+  func_shell5 = Parameter(2);
+  func_shell6 = Parameter(2);
+  func_shell7 = Parameter(2);
+  func_shell8 = Parameter(2);
+  func_shell9 = Parameter(2);
+  func_shell10 = Parameter(2);
 }
 
@@ -107,180 +335,180 @@
  */
 double OnionModel :: operator()(double q) {
-	double dp[56];
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = n_shells();
-	dp[1] = scale();
-	dp[2] = rad_core0();
-	dp[3] = sld_core0();
-	dp[4] = sld_solv();
-	dp[5] = 0.0;
-
-	dp[6] = sld_out_shell1();
-	dp[7] = sld_out_shell2();
-	dp[8] = sld_out_shell3();
-	dp[9] = sld_out_shell4();
-	dp[10] = sld_out_shell5();
-	dp[11] = sld_out_shell6();
-	dp[12] = sld_out_shell7();
-	dp[13] = sld_out_shell8();
-	dp[14] = sld_out_shell9();
-	dp[15] = sld_out_shell10();
-
-	dp[16] = sld_in_shell1();
-	dp[17] = sld_in_shell2();
-	dp[18] = sld_in_shell3();
-	dp[19] = sld_in_shell4();
-	dp[20] = sld_in_shell5();
-	dp[21] = sld_in_shell6();
-	dp[22] = sld_in_shell7();
-	dp[23] = sld_in_shell8();
-	dp[24] = sld_in_shell9();
-	dp[25] = sld_in_shell10();
-
-	dp[26] = A_shell1();
-	dp[27] = A_shell2();
-	dp[28] = A_shell3();
-	dp[29] = A_shell4();
-	dp[30] = A_shell5();
-	dp[31] = A_shell6();
-	dp[32] = A_shell7();
-	dp[33] = A_shell8();
-	dp[34] = A_shell9();
-	dp[35] = A_shell10();
-
-	dp[36] = thick_shell1();
-	dp[37] = thick_shell2();
-	dp[38] = thick_shell3();
-	dp[39] = thick_shell4();
-	dp[40] = thick_shell5();
-	dp[41] = thick_shell6();
-	dp[42] = thick_shell7();
-	dp[43] = thick_shell8();
-	dp[44] = thick_shell9();
-	dp[45] = thick_shell10();
-
-	dp[46] = func_shell1();
-	dp[47] = func_shell2();
-	dp[48] = func_shell3();
-	dp[49] = func_shell4();
-	dp[50] = func_shell5();
-	dp[51] = func_shell6();
-	dp[52] = func_shell7();
-	dp[53] = func_shell8();
-	dp[54] = func_shell9();
-	dp[55] = func_shell10();
-
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	rad_core0.get_weights(weights_rad);
-
-	// Get the dispersion points for the thick 1
-	vector<WeightPoint> weights_s1;
-	thick_shell1.get_weights(weights_s1);
-
-	// Get the dispersion points for the thick 2
-	vector<WeightPoint> weights_s2;
-	thick_shell2.get_weights(weights_s2);
-
-	// Get the dispersion points for the thick 3
-	vector<WeightPoint> weights_s3;
-	thick_shell3.get_weights(weights_s3);
-
-	// Get the dispersion points for the thick 4
-	vector<WeightPoint> weights_s4;
-	thick_shell4.get_weights(weights_s4);
-
-	// Get the dispersion points for the thick 5
-	vector<WeightPoint> weights_s5;
-	thick_shell5.get_weights(weights_s5);
-
-	// Get the dispersion points for the thick 6
-	vector<WeightPoint> weights_s6;
-	thick_shell6.get_weights(weights_s6);
-
-	// Get the dispersion points for the thick 7
-	vector<WeightPoint> weights_s7;
-	thick_shell7.get_weights(weights_s7);
-
-	// Get the dispersion points for the thick 8
-	vector<WeightPoint> weights_s8;
-	thick_shell8.get_weights(weights_s8);
-	// Get the dispersion points for the thick 9
-	vector<WeightPoint> weights_s9;
-	thick_shell9.get_weights(weights_s9);
-
-	// Get the dispersion points for the thick 10
-	vector<WeightPoint> weights_s10;
-	thick_shell10.get_weights(weights_s10);
-
-
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp[2] = weights_rad[i].value;
-		// Loop over radius weight points
-		for(size_t j=0; j<weights_s1.size(); j++) {
-			dp[36] = weights_s1[j].value;
-			// Loop over radius weight points
-			for(size_t k=0; k<weights_s2.size(); k++) {
-				dp[37] = weights_s2[k].value;
-				// Loop over radius weight points
-				for(size_t l=0; l<weights_s3.size(); l++) {
-					dp[38] = weights_s3[l].value;
-					// Loop over radius weight points
-					for(size_t m=0; m<weights_s4.size(); m++) {
-						dp[39] = weights_s4[m].value;
-						for(size_t n=0; n<weights_s5.size(); n++) {
-							dp[40] = weights_s5[n].value;
-							for(size_t o=0; o<weights_s6.size(); o++) {
-								dp[41] = weights_s6[o].value;
-								for(size_t p=0; p<weights_s7.size(); p++) {
-									dp[42] = weights_s7[p].value;
-									for(size_t t=0; t<weights_s8.size(); t++) {
-										dp[43] = weights_s8[t].value;
-										for(size_t r=0; r<weights_s9.size(); r++) {
-											dp[44] = weights_s9[r].value;
-											for(size_t s=0; s<weights_s10.size(); s++) {
-												dp[45] = weights_s10[s].value;
-												//Un-normalize Shells by volume
-												sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
-														*weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
-														*weights_s9[r].weight*weights_s10[s].weight
-														* so_kernel(dp,q) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
-																+weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
-																+weights_s9[r].value+weights_s10[s].value),3.0);
-												//Find average volume
-												vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
-													*weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
-													*weights_s9[r].weight*weights_s10[s].weight
-													* pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
-															+weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
-															+weights_s9[r].value+weights_s10[s].value),3.0);
-												norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
-													*weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
-													*weights_s9[r].weight*weights_s10[s].weight;
-											}
-										}
-									}
-								}
-							}
-						}
-					}
-				}
-			}
-		}
-	}
-
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[56];
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = n_shells();
+  dp[1] = scale();
+  dp[2] = rad_core0();
+  dp[3] = sld_core0();
+  dp[4] = sld_solv();
+  dp[5] = 0.0;
+
+  dp[6] = sld_out_shell1();
+  dp[7] = sld_out_shell2();
+  dp[8] = sld_out_shell3();
+  dp[9] = sld_out_shell4();
+  dp[10] = sld_out_shell5();
+  dp[11] = sld_out_shell6();
+  dp[12] = sld_out_shell7();
+  dp[13] = sld_out_shell8();
+  dp[14] = sld_out_shell9();
+  dp[15] = sld_out_shell10();
+
+  dp[16] = sld_in_shell1();
+  dp[17] = sld_in_shell2();
+  dp[18] = sld_in_shell3();
+  dp[19] = sld_in_shell4();
+  dp[20] = sld_in_shell5();
+  dp[21] = sld_in_shell6();
+  dp[22] = sld_in_shell7();
+  dp[23] = sld_in_shell8();
+  dp[24] = sld_in_shell9();
+  dp[25] = sld_in_shell10();
+
+  dp[26] = A_shell1();
+  dp[27] = A_shell2();
+  dp[28] = A_shell3();
+  dp[29] = A_shell4();
+  dp[30] = A_shell5();
+  dp[31] = A_shell6();
+  dp[32] = A_shell7();
+  dp[33] = A_shell8();
+  dp[34] = A_shell9();
+  dp[35] = A_shell10();
+
+  dp[36] = thick_shell1();
+  dp[37] = thick_shell2();
+  dp[38] = thick_shell3();
+  dp[39] = thick_shell4();
+  dp[40] = thick_shell5();
+  dp[41] = thick_shell6();
+  dp[42] = thick_shell7();
+  dp[43] = thick_shell8();
+  dp[44] = thick_shell9();
+  dp[45] = thick_shell10();
+
+  dp[46] = func_shell1();
+  dp[47] = func_shell2();
+  dp[48] = func_shell3();
+  dp[49] = func_shell4();
+  dp[50] = func_shell5();
+  dp[51] = func_shell6();
+  dp[52] = func_shell7();
+  dp[53] = func_shell8();
+  dp[54] = func_shell9();
+  dp[55] = func_shell10();
+
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  rad_core0.get_weights(weights_rad);
+
+  // Get the dispersion points for the thick 1
+  vector<WeightPoint> weights_s1;
+  thick_shell1.get_weights(weights_s1);
+
+  // Get the dispersion points for the thick 2
+  vector<WeightPoint> weights_s2;
+  thick_shell2.get_weights(weights_s2);
+
+  // Get the dispersion points for the thick 3
+  vector<WeightPoint> weights_s3;
+  thick_shell3.get_weights(weights_s3);
+
+  // Get the dispersion points for the thick 4
+  vector<WeightPoint> weights_s4;
+  thick_shell4.get_weights(weights_s4);
+
+  // Get the dispersion points for the thick 5
+  vector<WeightPoint> weights_s5;
+  thick_shell5.get_weights(weights_s5);
+
+  // Get the dispersion points for the thick 6
+  vector<WeightPoint> weights_s6;
+  thick_shell6.get_weights(weights_s6);
+
+  // Get the dispersion points for the thick 7
+  vector<WeightPoint> weights_s7;
+  thick_shell7.get_weights(weights_s7);
+
+  // Get the dispersion points for the thick 8
+  vector<WeightPoint> weights_s8;
+  thick_shell8.get_weights(weights_s8);
+  // Get the dispersion points for the thick 9
+  vector<WeightPoint> weights_s9;
+  thick_shell9.get_weights(weights_s9);
+
+  // Get the dispersion points for the thick 10
+  vector<WeightPoint> weights_s10;
+  thick_shell10.get_weights(weights_s10);
+
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[2] = weights_rad[i].value;
+    // Loop over radius weight points
+    for(size_t j=0; j<weights_s1.size(); j++) {
+      dp[36] = weights_s1[j].value;
+      // Loop over radius weight points
+      for(size_t k=0; k<weights_s2.size(); k++) {
+        dp[37] = weights_s2[k].value;
+        // Loop over radius weight points
+        for(size_t l=0; l<weights_s3.size(); l++) {
+          dp[38] = weights_s3[l].value;
+          // Loop over radius weight points
+          for(size_t m=0; m<weights_s4.size(); m++) {
+            dp[39] = weights_s4[m].value;
+            for(size_t n=0; n<weights_s5.size(); n++) {
+              dp[40] = weights_s5[n].value;
+              for(size_t o=0; o<weights_s6.size(); o++) {
+                dp[41] = weights_s6[o].value;
+                for(size_t p=0; p<weights_s7.size(); p++) {
+                  dp[42] = weights_s7[p].value;
+                  for(size_t t=0; t<weights_s8.size(); t++) {
+                    dp[43] = weights_s8[t].value;
+                    for(size_t r=0; r<weights_s9.size(); r++) {
+                      dp[44] = weights_s9[r].value;
+                      for(size_t s=0; s<weights_s10.size(); s++) {
+                        dp[45] = weights_s10[s].value;
+                        //Un-normalize Shells by volume
+                        sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
+                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
+                            *weights_s9[r].weight*weights_s10[s].weight
+                            * so_kernel(dp,q) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
+                                +weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
+                                +weights_s9[r].value+weights_s10[s].value),3.0);
+                        //Find average volume
+                        vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
+                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
+                            *weights_s9[r].weight*weights_s10[s].weight
+                            * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
+                                +weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
+                                +weights_s9[r].value+weights_s10[s].value),3.0);
+                        norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
+                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
+                            *weights_s9[r].weight*weights_s10[s].weight;
+                      }
+                    }
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -292,6 +520,6 @@
  */
 double OnionModel :: operator()(double qx, double qy) {
-	double q = sqrt(qx*qx + qy*qy);
-	return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -304,5 +532,5 @@
  */
 double OnionModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -312,127 +540,126 @@
  */
 double OnionModel :: calculate_ER() {
-	OnionParameters dp;
-	dp.rad_core0 = rad_core0();
-	dp.thick_shell1 = thick_shell1();
-	dp.thick_shell2 = thick_shell2();
-	dp.thick_shell3 = thick_shell3();
-	dp.thick_shell4 = thick_shell4();
-	dp.thick_shell5 = thick_shell5();
-	dp.thick_shell6 = thick_shell6();
-	dp.thick_shell7 = thick_shell7();
-	dp.thick_shell8 = thick_shell8();
-	dp.thick_shell9 = thick_shell9();
-	dp.thick_shell10 = thick_shell10();
-
-
-	double rad_out = 0.0;
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	rad_core0.get_weights(weights_rad);
-
-	// Get the dispersion points for the thick 1
-	vector<WeightPoint> weights_s1;
-	thick_shell1.get_weights(weights_s1);
-
-	// Get the dispersion points for the thick 2
-	vector<WeightPoint> weights_s2;
-	thick_shell2.get_weights(weights_s2);
-
-	// Get the dispersion points for the thick 3
-	vector<WeightPoint> weights_s3;
-	thick_shell3.get_weights(weights_s3);
-
-	// Get the dispersion points for the thick 4
-	vector<WeightPoint> weights_s4;
-	thick_shell4.get_weights(weights_s4);
-	// Get the dispersion points for the thick 5
-	vector<WeightPoint> weights_s5;
-	thick_shell5.get_weights(weights_s5);
-
-	// Get the dispersion points for the thick 6
-	vector<WeightPoint> weights_s6;
-	thick_shell6.get_weights(weights_s6);
-
-	// Get the dispersion points for the thick 7
-	vector<WeightPoint> weights_s7;
-	thick_shell7.get_weights(weights_s7);
-
-	// Get the dispersion points for the thick 8
-	vector<WeightPoint> weights_s8;
-	thick_shell8.get_weights(weights_s8);
-	// Get the dispersion points for the thick 9
-	vector<WeightPoint> weights_s9;
-	thick_shell9.get_weights(weights_s9);
-
-	// Get the dispersion points for the thick 10
-	vector<WeightPoint> weights_s10;
-	thick_shell10.get_weights(weights_s10);
-
-
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp.rad_core0 = weights_rad[i].value;
-		// Loop over radius weight points
-		for(size_t j=0; j<weights_s1.size(); j++) {
-			dp.thick_shell1 = weights_s1[j].value;
-			// Loop over radius weight points
-			for(size_t k=0; k<weights_s2.size(); k++) {
-				dp.thick_shell2 = weights_s2[k].value;
-				// Loop over radius weight points
-				for(size_t l=0; l<weights_s3.size(); l++) {
-					dp.thick_shell3 = weights_s3[l].value;
-					// Loop over radius weight points
-					for(size_t m=0; m<weights_s4.size(); m++) {
-						dp.thick_shell4 = weights_s4[m].value;
-						// Loop over radius weight points
-						for(size_t n=0; j<weights_s5.size(); n++) {
-							dp.thick_shell5 = weights_s5[n].value;
-							// Loop over radius weight points
-							for(size_t o=0; k<weights_s6.size(); o++) {
-								dp.thick_shell6 = weights_s6[o].value;
-								// Loop over radius weight points
-								for(size_t p=0; l<weights_s7.size(); p++) {
-									dp.thick_shell7 = weights_s7[p].value;
-									// Loop over radius weight points
-									for(size_t t=0; m<weights_s8.size(); t++) {
-										dp.thick_shell8 = weights_s8[t].value;
-										// Loop over radius weight points
-										for(size_t r=0; l<weights_s9.size(); r++) {
-											dp.thick_shell8 = weights_s9[r].value;
-											// Loop over radius weight points
-											for(size_t s=0; m<weights_s10.size(); s++) {
-												dp.thick_shell10 = weights_s10[s].value;
-												//Un-normalize FourShell by volume
-												sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
-													*weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
-													*weights_s9[r].weight*weights_s10[s].weight
-													* (dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4+dp.thick_shell5
-															+dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10);
-												norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight
-													*weights_s4[m].weight*weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight
-													*weights_s8[t].weight*weights_s9[r].weight*weights_s10[s].weight;
-											}
-										}
-									}
-								}
-							}
-						}
-					}
-				}
-			}
-		}
-	}
-
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		rad_out = dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4
-					+dp.thick_shell5+dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10;}
-	return rad_out;
+  OnionParameters dp;
+  dp.rad_core0 = rad_core0();
+  dp.thick_shell1 = thick_shell1();
+  dp.thick_shell2 = thick_shell2();
+  dp.thick_shell3 = thick_shell3();
+  dp.thick_shell4 = thick_shell4();
+  dp.thick_shell5 = thick_shell5();
+  dp.thick_shell6 = thick_shell6();
+  dp.thick_shell7 = thick_shell7();
+  dp.thick_shell8 = thick_shell8();
+  dp.thick_shell9 = thick_shell9();
+  dp.thick_shell10 = thick_shell10();
+
+  double rad_out = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  rad_core0.get_weights(weights_rad);
+
+  // Get the dispersion points for the thick 1
+  vector<WeightPoint> weights_s1;
+  thick_shell1.get_weights(weights_s1);
+
+  // Get the dispersion points for the thick 2
+  vector<WeightPoint> weights_s2;
+  thick_shell2.get_weights(weights_s2);
+
+  // Get the dispersion points for the thick 3
+  vector<WeightPoint> weights_s3;
+  thick_shell3.get_weights(weights_s3);
+
+  // Get the dispersion points for the thick 4
+  vector<WeightPoint> weights_s4;
+  thick_shell4.get_weights(weights_s4);
+  // Get the dispersion points for the thick 5
+  vector<WeightPoint> weights_s5;
+  thick_shell5.get_weights(weights_s5);
+
+  // Get the dispersion points for the thick 6
+  vector<WeightPoint> weights_s6;
+  thick_shell6.get_weights(weights_s6);
+
+  // Get the dispersion points for the thick 7
+  vector<WeightPoint> weights_s7;
+  thick_shell7.get_weights(weights_s7);
+
+  // Get the dispersion points for the thick 8
+  vector<WeightPoint> weights_s8;
+  thick_shell8.get_weights(weights_s8);
+  // Get the dispersion points for the thick 9
+  vector<WeightPoint> weights_s9;
+  thick_shell9.get_weights(weights_s9);
+
+  // Get the dispersion points for the thick 10
+  vector<WeightPoint> weights_s10;
+  thick_shell10.get_weights(weights_s10);
+
+
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp.rad_core0 = weights_rad[i].value;
+    // Loop over radius weight points
+    for(size_t j=0; j<weights_s1.size(); j++) {
+      dp.thick_shell1 = weights_s1[j].value;
+      // Loop over radius weight points
+      for(size_t k=0; k<weights_s2.size(); k++) {
+        dp.thick_shell2 = weights_s2[k].value;
+        // Loop over radius weight points
+        for(size_t l=0; l<weights_s3.size(); l++) {
+          dp.thick_shell3 = weights_s3[l].value;
+          // Loop over radius weight points
+          for(size_t m=0; m<weights_s4.size(); m++) {
+            dp.thick_shell4 = weights_s4[m].value;
+            // Loop over radius weight points
+            for(size_t n=0; j<weights_s5.size(); n++) {
+              dp.thick_shell5 = weights_s5[n].value;
+              // Loop over radius weight points
+              for(size_t o=0; k<weights_s6.size(); o++) {
+                dp.thick_shell6 = weights_s6[o].value;
+                // Loop over radius weight points
+                for(size_t p=0; l<weights_s7.size(); p++) {
+                  dp.thick_shell7 = weights_s7[p].value;
+                  // Loop over radius weight points
+                  for(size_t t=0; m<weights_s8.size(); t++) {
+                    dp.thick_shell8 = weights_s8[t].value;
+                    // Loop over radius weight points
+                    for(size_t r=0; l<weights_s9.size(); r++) {
+                      dp.thick_shell8 = weights_s9[r].value;
+                      // Loop over radius weight points
+                      for(size_t s=0; m<weights_s10.size(); s++) {
+                        dp.thick_shell10 = weights_s10[s].value;
+                        //Un-normalize FourShell by volume
+                        sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
+                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
+                            *weights_s9[r].weight*weights_s10[s].weight
+                            * (dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4+dp.thick_shell5
+                                +dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10);
+                        norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight
+                            *weights_s4[m].weight*weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight
+                            *weights_s8[t].weight*weights_s9[r].weight*weights_s10[s].weight;
+                      }
+                    }
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    rad_out = dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4
+        +dp.thick_shell5+dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10;}
+  return rad_out;
 }
Index: sansmodels/src/python_wrapper/COnionModel.cpp
===================================================================
--- sansmodels/src/python_wrapper/COnionModel.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/python_wrapper/COnionModel.cpp	(revision c637521b288a6e3976fc9c0d461b44f70bb58c10)
@@ -18,5 +18,5 @@
  *
  * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
- *          DO NOT MODIFY THIS FILE, MODIFY onion.h
+ *          DO NOT MODIFY THIS FILE, MODIFY ../c_extensions/onion.h
  *          AND RE-RUN THE GENERATOR SCRIPT
  *
@@ -33,8 +33,8 @@
 #include <math.h>
 #include <time.h>
+
+}
+
 #include "onion.h"
-}
-
-#include "models.hh"
 #include "dispersion_visitor.hh"
 
