Index: sansmodels/src/c_extensions/DiamCyl.h
===================================================================
--- sansmodels/src/c_extensions/DiamCyl.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/DiamCyl.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,4 +1,5 @@
 #if !defined(DiamCyl_h)
 #define DiamCyl_h
+#include "parameters.hh"
 
 /**
@@ -16,23 +17,23 @@
 **/
 
-typedef struct {
-    ///	Radius [A]
-    //  [DEFAULT]=radius=20.0 A
-    double radius;
+class DiamCylFunc{
+public:
+  // Model parameters
+  /// Radius [A]
+  //  [DEFAULT]=radius=20.0 A
+  Parameter radius;
+  /// Length [A]
+  //  [DEFAULT]=length= 400 A
+  Parameter length;
 
-    /// Length [A]
-    //  [DEFAULT]=length= 400 A
-    double length;
+  // Constructor
+  DiamCylFunc();
 
-} DiamCyldParameters;
-
-
-
-/// 1D scattering function
-//double DiamCyld_analytical_1D(DiamCyldParameters *pars, double q);
-
-/// 2D scattering function
-//double DiamCyld_analytical_2D(DiamCyldParameters *pars, double q, double phi);
-//double DiamCyld_analytical_2DXY(DiamCyldParameters *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_extensions/DiamEllip.h
===================================================================
--- sansmodels/src/c_extensions/DiamEllip.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/DiamEllip.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,4 +1,5 @@
 #if !defined(DiamEllip_h)
 #define DiamEllip_h
+#include "parameters.hh"
 
 /**
@@ -20,23 +21,25 @@
 
  **/
-typedef struct {
-    ///	Polar radius [A]
-    //  [DEFAULT]=radius_a=20.0 A
-    double radius_a;
 
-    /// Equatorial radius [A]
-    //  [DEFAULT]=radius_b= 400 A
-    double radius_b;
+class DiamEllipFunc{
+public:
+  // Model parameters
+  /// Polar radius [A]
+  //  [DEFAULT]=radius_a=20.0 A
+  Parameter radius_a;
 
-} DiamEllipsParameters;
+  /// Equatorial radius [A]
+  //  [DEFAULT]=radius_b= 400 A
+  Parameter radius_b;
 
+  // Constructor
+  DiamEllipFunc();
 
-
-/// 1D scattering function
-//double DiamEllips_analytical_1D(DiamEllipsParameters *pars, double q);
-
-/// 2D scattering function
-//double DiamEllips_analytical_2D(DiamEllipsParameters *pars, double q, double phi);
-//double DiamEllips_analytical_2DXY(DiamEllipsParameters *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_extensions/Hardsphere.h
===================================================================
--- sansmodels/src/c_extensions/Hardsphere.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/Hardsphere.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,44 +1,46 @@
 #if !defined(Hardsphere_h)
 #define Hardsphere_h
+#include "parameters.hh"
 
 /**
  * Structure definition for HardsphereStructure (factor) parameters
  */
- //[PYTHONCLASS] = HardsphereStructure
- //[DISP_PARAMS] = effect_radius
- //[DESCRIPTION] =<text>Structure factor for interacting particles:                   .
- //
- //  The interparticle potential is
- //
- //                     U(r)= inf   , r < 2R
- //                         = 0     , r >= 2R
- //
- //						R: effective radius of the Hardsphere particle
- //						V:The volume fraction
- //
- //    Ref: Percus., J. K.,etc., J. Phy.
- //     Rev. 1958, 110, 1.
- //	 </text>
- //[FIXED]= effect_radius.width
+//[PYTHONCLASS] = HardsphereStructure
+//[DISP_PARAMS] = effect_radius
+//[DESCRIPTION] =<text>Structure factor for interacting particles:                   .
+//
+//  The interparticle potential is
+//
+//                     U(r)= inf   , r < 2R
+//                         = 0     , r >= 2R
+//
+//						R: effective radius of the Hardsphere particle
+//						V:The volume fraction
+//
+//    Ref: Percus., J. K.,etc., J. Phy.
+//     Rev. 1958, 110, 1.
+//	 </text>
+//[FIXED]= effect_radius.width
 
+class HardsphereStructure{
+public:
+  // Model parameters
+  /// effect radius of hardsphere [A]
+  //  [DEFAULT]=effect_radius=50.0 [A]
+  Parameter effect_radius;
 
-typedef struct {
-    /// effect radius of hardsphere [A]
-    //  [DEFAULT]=effect_radius=50.0 [A]
-    double effect_radius;
+  /// Volume fraction
+  //  [DEFAULT]=volfraction= 0.2
+  Parameter volfraction;
 
-    ///	Volume fraction
-    //  [DEFAULT]=volfraction= 0.2
-    double volfraction;
-} HardsphereParameters;
+  // Constructor
+  HardsphereStructure();
 
-
-
-/// 1D scattering function
-//double Hardsphere_analytical_1D(HardsphereParameters *pars, double q);
-
-/// 2D scattering function
-//double Hardsphere_analytical_2D(HardsphereParameters *pars, double q, double phi);
-//double Hardsphere_analytical_2DXY(HardsphereParameters *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_extensions/HayterMSA.h
===================================================================
--- sansmodels/src/c_extensions/HayterMSA.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/HayterMSA.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,59 +1,64 @@
 #if !defined(HayterMSA_h)
 #define HayterMSA_h
+#include "parameters.hh"
 
 /**
  * Structure definition for screened Coulomb interaction
  */
- //[PYTHONCLASS] = HayterMSAStructure
- //[DISP_PARAMS] = effect_radius
- //[DESCRIPTION] =<text>To calculate the structure factor (the Fourier transform of the
- //                     pair correlation function g(r)) for a system of
- //                     charged, spheroidal objects in a dielectric
- //						medium.
- //                     When combined with an appropriate form
- //                     factor, this allows for inclusion of
- //                     the interparticle interference effects
- //                     due to screened coulomb repulsion between
- //                     charged particles.
- //						(Note: charge > 0 required.)
- //
- //                     Ref: JP Hansen and JB Hayter, Molecular
- //                           Physics 46, 651-656 (1982).
- //
- //				</text>
- //[FIXED]= effect_radius.width
+//[PYTHONCLASS] = HayterMSAStructure
+//[DISP_PARAMS] = effect_radius
+//[DESCRIPTION] =<text>To calculate the structure factor (the Fourier transform of the
+//                     pair correlation function g(r)) for a system of
+//                     charged, spheroidal objects in a dielectric
+//						medium.
+//                     When combined with an appropriate form
+//                     factor, this allows for inclusion of
+//                     the interparticle interference effects
+//                     due to screened coulomb repulsion between
+//                     charged particles.
+//						(Note: charge > 0 required.)
+//
+//                     Ref: JP Hansen and JB Hayter, Molecular
+//                           Physics 46, 651-656 (1982).
+//
+//				</text>
+//[FIXED]= effect_radius.width
 
-typedef struct {
-    ///	effetitve radius of particle [A]
-    //  [DEFAULT]=effect_radius=20.75 [A]
-    double effect_radius;
+class HayterMSAStructure{
+public:
+  // Model parameters
+  /// effetitve radius of particle [A]
+  //  [DEFAULT]=effect_radius=20.75 [A]
+  Parameter effect_radius;
 
-    /// charge
-    //  [DEFAULT]=charge= 19
-    double charge;
+  /// charge
+  //  [DEFAULT]=charge= 19
+  Parameter charge;
 
-    /// Volume fraction
-    //  [DEFAULT]=volfraction= 0.0192
-    double volfraction;
+  /// Volume fraction
+  //  [DEFAULT]=volfraction= 0.0192
+  Parameter volfraction;
 
-	///	Temperature [K]
-    //  [DEFAULT]=temperature= 318.16 [K]
-    double temperature;
+  /// Temperature [K]
+  //  [DEFAULT]=temperature= 318.16 [K]
+  Parameter temperature;
 
-	///	Monovalent salt concentration [M]
-    //  [DEFAULT]=saltconc= 0 [M]
-    double saltconc;
+  /// Monovalent salt concentration [M]
+  //  [DEFAULT]=saltconc= 0 [M]
+  Parameter saltconc;
 
-    ///	Dielectric constant of solvent
-    //  [DEFAULT]=dielectconst= 71.08
-    double dielectconst;
-} HayterMSAParameters;
+  /// Dielectric constant of solvent
+  //  [DEFAULT]=dielectconst= 71.08
+  Parameter dielectconst;
 
-/// 1D scattering function
-//double HayterMSA_analytical_1D(HayterMSAParameters *pars, double q);
+  // Constructor
+  HayterMSAStructure();
 
-/// 2D scattering function
-//double HayterMSA_analytical_2D(HayterMSAParameters *pars, double q, double phi);
-//double HayterMSA_analytical_2DXY(HayterMSAParameters *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_extensions/SquareWell.h
===================================================================
--- sansmodels/src/c_extensions/SquareWell.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/SquareWell.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,7 +1,8 @@
 #if !defined(SquareWell_h)
 #define SquareWell_h
+#include "parameters.hh"
 
 /**Structure definition for SquareWell parameters
-*/
+ */
 //   [PYTHONCLASS] = SquareWellStructure
 //   [DISP_PARAMS] = effect_radius
@@ -26,32 +27,32 @@
 //[ORIENTATION_PARAMS]= <text> </text>
 
+class SquareWellStructure{
+public:
+  // Model parameters
+  /// effective radius of particle [A]
+  //  [DEFAULT]=effect_radius=50.0 [A]
+  Parameter effect_radius;
 
-typedef struct {
-    ///	effective radius of particle [A]
-    //  [DEFAULT]=effect_radius=50.0 [A]
-    double effect_radius;
+  /// Volume fraction
+  //  [DEFAULT]=volfraction= 0.040
+  Parameter volfraction;
 
-    /// Volume fraction
-    //  [DEFAULT]=volfraction= 0.040
-    double volfraction;
+  /// Well depth [kT]
+  //  [DEFAULT]=welldepth= 1.50 [kT]
+  Parameter welldepth;
 
-    ///	Well depth [kT]
-    //  [DEFAULT]=welldepth= 1.50 [kT]
-    double welldepth;
+  /// Well width
+  //  [DEFAULT]=wellwidth= 1.20
+  Parameter wellwidth;
 
-    ///	Well width
-    //  [DEFAULT]=wellwidth= 1.20
-    double wellwidth;
+  // Constructor
+  SquareWellStructure();
 
-} SquareWellParameters;
-
-
-
-/// 1D scattering function
-//double SquareWell_analytical_1D(SquareWellParameters *pars, double q);
-
-/// 2D scattering function
-//double SquareWell_analytical_2D(SquareWellParameters *pars, double q, double phi);
-//double SquareWell_analytical_2DXY(SquareWellParameters *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_extensions/StickyHS.h
===================================================================
--- sansmodels/src/c_extensions/StickyHS.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/StickyHS.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,54 +1,58 @@
 #if !defined(StickyHS_h)
 #define StickyHS_h
+#include "parameters.hh"
 
 /**
  * Structure definition for StickyHS_struct (struncture factor) parameters
  */
- //[PYTHONCLASS] = StickyHSStructure
- //[DISP_PARAMS] = effect_radius
- //[DESCRIPTION] =<text> Structure Factor for interacting particles:                               .
- //
- //  The interaction potential is
- //
- //                     U(r)= inf , r < 2R
- //                         = -Uo  , 2R < r < 2R + w
- //                         = 0   , r >= 2R +w
- //
- //						R: effective radius of the hardsphere
- //                     stickiness = [exp(Uo/kT)]/(12*perturb)
- //                     perturb = w/(w+ 2R) , 0.01 =< w <= 0.1
- //                     w: The width of the square well ,w > 0
- //						v: The volume fraction , v > 0
- //
- //                     Ref: Menon, S. V. G.,et.al., J. Chem.
- //                      Phys., 1991, 95(12), 9186-9190.
- //				</text>
- //[FIXED]= effect_radius.width
-typedef struct {
-    /// effect radius of hardsphere [A]
-    //  [DEFAULT]=effect_radius=50.0 [A]
-    double effect_radius;
+//[PYTHONCLASS] = StickyHSStructure
+//[DISP_PARAMS] = effect_radius
+//[DESCRIPTION] =<text> Structure Factor for interacting particles:                               .
+//
+//  The interaction potential is
+//
+//                     U(r)= inf , r < 2R
+//                         = -Uo  , 2R < r < 2R + w
+//                         = 0   , r >= 2R +w
+//
+//						R: effective radius of the hardsphere
+//                     stickiness = [exp(Uo/kT)]/(12*perturb)
+//                     perturb = w/(w+ 2R) , 0.01 =< w <= 0.1
+//                     w: The width of the square well ,w > 0
+//						v: The volume fraction , v > 0
+//
+//                     Ref: Menon, S. V. G.,et.al., J. Chem.
+//                      Phys., 1991, 95(12), 9186-9190.
+//				</text>
+//[FIXED]= effect_radius.width
 
-    ///	Volume fraction
-    //  [DEFAULT]=volfraction= 0.1
-    double volfraction;
+class StickyHSStructure{
+public:
+  // Model parameters
+  /// effect radius of hardsphere [A]
+  //  [DEFAULT]=effect_radius=50.0 [A]
+  Parameter effect_radius;
 
-	/// Perturbation Parameter ( between 0.01 - 0.1)
-    //  [DEFAULT]=perturb=0.05
-    double perturb;
+  /// Volume fraction
+  //  [DEFAULT]=volfraction= 0.1
+  Parameter volfraction;
 
-    ///	Stickiness
-    //  [DEFAULT]=stickiness= 0.2
-    double stickiness;
-} StickyHSParameters;
+  /// Perturbation Parameter ( between 0.01 - 0.1)
+  //  [DEFAULT]=perturb=0.05
+  Parameter perturb;
 
+  /// Stickiness
+  //  [DEFAULT]=stickiness= 0.2
+  Parameter stickiness;
 
+  // Constructor
+  StickyHSStructure();
 
-/// 1D scattering function
-//double StickyHS_analytical_1D(StickyHSParameters *pars, double q);
-
-/// 2D scattering function
-//double StickyHS_analytical_2D(StickyHSParameters *pars, double q, double phi);
-//double StickyHS_analytical_2DXY(StickyHSParameters *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_extensions/bcc.c
===================================================================
--- sansmodels/src/c_extensions/bcc.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,146 +1,0 @@
-/*
- * Scattering model for a BC_ParaCrystal
- */
-#include "bcc.h"
-#include "libSphere.h"
-#include <math.h>
-#include <stdio.h>
-
-
-/**
- * Function to evaluate 1D scattering function
- * @param pars: parameters of the BCC_ParaCrystal
- * @param q: q-value
- * @return: function value
- */
-double bcc_analytical_1D(BCParameters *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 = BCC_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 BCC_ParaCrystal
- * @param q: q-value
- * @return: function value
- */
-double bc_analytical_2DXY(BCParameters *pars, double qx, double qy){
-	double q;
-	q = sqrt(qx*qx+qy*qy);
-	return bc_analytical_2D_scaled(pars, q, qx/q, qy/q);
-}
-
-double bc_analytical_2D(BCParameters *pars, double q, double phi) {
-	return bc_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
-}
-
-/**
- * Function to evaluate 2D scattering function
- * @param pars: parameters of the BCCCrystalModel
- * @param q: q-value
- * @param q_x: q_x / q
- * @param q_y: q_y / q
- * @return: function value
- */
-double bc_analytical_2D_scaled(BCParameters *pars, double q, double q_x, double q_y) {
-	double b3_x, b3_y, b3_z, b1_x, b1_y;
-	double q_z;
-	double alpha, cos_val_b3, cos_val_b2, cos_val_b1;
-	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, Fkq, Fkq_2;
-
-  //convert angle degree to radian
-  double pi = 4.0*atan(1.0);
-  double theta = pars->theta * pi/180.0;
-  double phi = pars->phi * pi/180.0;
-  double psi = pars->psi * pi/180.0;
-
-  double dp[5];
-	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);
-
-	//the occupied volume of the lattice
-	latticeScale = 2.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa/sqrt(3.0/4.0),3.0);
-	// q vector
-	q_z = 0.0; // for SANS; assuming qz is negligible
-	/// Angles here are respect to detector coordinate
-	///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
-    // b3 axis orientation
-    b3_x = sin(theta) * cos(phi);//negative sign here???
-    b3_y = sin(theta) * sin(phi);
-    b3_z = cos(theta);
-    cos_val_b3 =  b3_x*q_x + b3_y*q_y + b3_z*q_z;
-
-    alpha = acos(cos_val_b3);
-    // b1 axis orientation
-    b1_x = sin(psi);
-    b1_y = cos(psi);
-	cos_val_b1 = (b1_x*q_x + b1_y*q_y);
-    // b2 axis orientation
-	cos_val_b2 = sin(acos(cos_val_b1));
-	// alpha corrections
-	cos_val_b2 *= sin(alpha);
-	cos_val_b1 *= sin(alpha);
-
-    // Compute the angle btw vector q and the a3 axis
-    a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1-cos_val_b3);
-
-    // a1 axis
-    a1_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b2-cos_val_b1);
-
-    // a2 axis
-    a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1-cos_val_b2);
-
-    // The following test should always pass
-    if (fabs(cos_val_b3)>1.0) {
-    	printf("bcc_ana_2D: Unexpected error: cos(alpha)>1\n");
-     	return 0;
-    }
-    // Get Fkq and Fkq_2
-    Fkq = exp(-0.5*pow(Da/aa,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
-    Fkq_2 = Fkq*Fkq;
-    // Call Zq=Z1*Z2*Z3
-    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
-    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_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/bcc.h
===================================================================
--- sansmodels/src/c_extensions/bcc.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/bcc.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,81 +1,83 @@
 #if !defined(bcc_h)
 #define bcc_h
+#include "parameters.hh"
 
 /**
  * Structure definition for BCC_ParaCrystal parameters
  */
- //[PYTHONCLASS] = BCCrystalModel
- //[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 face centered cubic structure.
- //				[Body Centered 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>
+//[PYTHONCLASS] = BCCrystalModel
+//[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 face centered cubic structure.
+//				[Body Centered 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;
+class BCCrystalModel{
+public:
+  // Model parameters
+  /// Scale factor
+  //  [DEFAULT]=scale= 1.0
+  Parameter scale;
 
-    ///	Nearest neighbor distance [A]
-    //  [DEFAULT]=dnn=220.0 [A]
-    double dnn;
+  /// Nearest neighbor distance [A]
+  //  [DEFAULT]=dnn=220.0 [A]
+  Parameter dnn;
 
-    ///	Paracrystal distortion factor g
-    //  [DEFAULT]=d_factor=0.06
-    double d_factor;
+  /// Paracrystal distortion factor g
+  //  [DEFAULT]=d_factor=0.06
+  Parameter d_factor;
 
-    ///	Radius of sphere [A]
-    //  [DEFAULT]=radius=40.0 [A]
-    double radius;
+  /// Radius of sphere [A]
+  //  [DEFAULT]=radius=40.0 [A]
+  Parameter radius;
 
-    ///	sldSph [1/A^(2)]
-    //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)]
-    double sldSph;
+  /// 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)]
-    double sldSolv;
+  /// sldSolv [1/A^(2)]
+  //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)]
+  Parameter 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;
+  /// 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;
 
-} BCParameters;
+  // Constructor
+  BCCrystalModel();
+
+  // 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);
+};
 
 
-
-/// 1D scattering function
-double bc_analytical_1D(BCParameters *pars, double q);
-
-/// 2D scattering function
-double bc_analytical_2D(BCParameters *pars, double q, double phi);
-double bc_analytical_2DXY(BCParameters *pars, double qx, double qy);
-double bc_analytical_2D_scaled(BCParameters *pars, double q, double q_x, double q_y);
-
 #endif
Index: sansmodels/src/c_extensions/fcc.c
===================================================================
--- sansmodels/src/c_extensions/fcc.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,142 +1,0 @@
-/*
- * Scattering model for a SC_ParaCrystal
- */
-#include "fcc.h"
-#include "libSphere.h"
-#include <math.h>
-#include <stdio.h>
-
-/**
- * Function to evaluate 1D scattering function
- * @param pars: parameters of the FCC_ParaCrystal
- * @param q: q-value
- * @return: function value
- */
-double fcc_analytical_1D(FCParameters *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 = FCC_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 FCC_ParaCrystal
- * @param q: q-value
- * @return: function value
- */
-double fc_analytical_2DXY(FCParameters *pars, double qx, double qy){
-	double q;
-	q = sqrt(qx*qx+qy*qy);
-	return fc_analytical_2D_scaled(pars, q, qx/q, qy/q);
-}
-
-double fc_analytical_2D(FCParameters *pars, double q, double phi) {
-	return fc_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
-}
-
-/**
- * Function to evaluate 2D scattering function
- * @param pars: parameters of the FCCCrystalModel
- * @param q: q-value
- * @param q_x: q_x / q
- * @param q_y: q_y / q
- * @return: function value
- */
-double fc_analytical_2D_scaled(FCParameters *pars, double q, double q_x, double q_y) {
-	double b3_x, b3_y, b3_z, b1_x, b1_y;
-	double q_z;
-	double alpha, cos_val_b3, cos_val_b2, cos_val_b1;
-	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, Fkq, Fkq_2;
-
-	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);
-	//contrast = pars->sldSph - pars->sldSolv;
-
-	latticeScale = 4.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa*sqrt(2.0),3.0);
-	// q vector
-	q_z = 0.0; // for SANS; assuming qz is negligible
-	/// Angles here are respect to detector coordinate
-	///  instead of against q coordinate in PRB 36(46), 3(6), 1754(3854)
-    // b3 axis orientation
-    b3_x = sin(theta) * cos(phi);//negative sign here???
-    b3_y = sin(theta) * sin(phi);
-    b3_z = cos(theta);
-    cos_val_b3 =  b3_x*q_x + b3_y*q_y + b3_z*q_z;
-    alpha = acos(cos_val_b3);
-    // b1 axis orientation
-    b1_x = sin(psi);
-    b1_y = cos(psi);
-	cos_val_b1 = (b1_x*q_x + b1_y*q_y);
-    // b2 axis orientation
-	cos_val_b2 = sin(acos(cos_val_b1));
-	// alpha correction
-	cos_val_b2 *= sin(alpha);
-	cos_val_b1 *= sin(alpha);
-
-    // Compute the angle btw vector q and the a3 axis
-    a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1);
-
-    // a1 axis
-    a1_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b3);
-
-    // a2 axis
-    a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1);
-
-    // The following test should always pass
-    if (fabs(cos_val_b3)>1.0) {
-    	printf("fcc_ana_2D: Unexpected error: cos(alpha)>1\n");
-     	return 0;
-    }
-    // Get Fkq and Fkq_2
-    Fkq = exp(-0.5*pow(Da/aa,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
-    Fkq_2 = Fkq*Fkq;
-    // Call Zq=Z1*Z2*Z3
-    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
-    Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
-    Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_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/fcc.h
===================================================================
--- sansmodels/src/c_extensions/fcc.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/fcc.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,4 +1,5 @@
 #if !defined(fcc_h)
 #define fcc_h
+#include "parameters.hh"
 
 /**
@@ -29,53 +30,53 @@
  //[ORIENTATION_PARAMS]= <text> phi;psi; theta; phi.width;psi.width; theta.width</text>
 
-typedef struct {
-    /// Scale factor
-    //  [DEFAULT]=scale= 1.0
-    double scale;
+class FCCrystalModel{
+public:
+  // Model parameters
+  /// Scale factor
+  //  [DEFAULT]=scale= 1.0
+  Parameter scale;
 
-    ///	Nearest neighbor distance [A]
-    //  [DEFAULT]=dnn=220.0 [A]
-    double dnn;
+  /// Nearest neighbor distance [A]
+  //  [DEFAULT]=dnn=220.0 [A]
+  Parameter dnn;
 
-    ///	Paracrystal distortion factor g
-    //  [DEFAULT]=d_factor=0.06
-    double d_factor;
+  /// Paracrystal distortion factor g
+  //  [DEFAULT]=d_factor=0.06
+  Parameter d_factor;
 
-    ///	Radius of sphere [A]
-    //  [DEFAULT]=radius=40.0 [A]
-    double radius;
+  /// Radius of sphere [A]
+  //  [DEFAULT]=radius=40.0 [A]
+  Parameter radius;
 
-    ///	sldSph [1/A^(2)]
-    //  [DEFAULT]=sldSph= 3.0e-6 [1/A^(2)]
-    double sldSph;
+  /// 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)]
-    double sldSolv;
+  /// sldSolv [1/A^(2)]
+  //  [DEFAULT]=sldSolv= 6.3e-6 [1/A^(2)]
+  Parameter 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;
+/// 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;
 
-} FCParameters;
+  // Constructor
+  FCCrystalModel();
 
-
-
-/// 1D scattering function
-double fc_analytical_1D(FCParameters *pars, double q);
-
-/// 2D scattering function
-double fc_analytical_2D(FCParameters *pars, double q, double phi);
-double fc_analytical_2DXY(FCParameters *pars, double qx, double qy);
-double fc_analytical_2D_scaled(FCParameters *pars, double q, double q_x, double q_y);
+  // 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_extensions/fuzzysphere.c
===================================================================
--- sansmodels/src/c_extensions/fuzzysphere.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,78 +1,0 @@
-/**
- * Scattering model for a fuzzy sphere
- */
-
-#include <math.h>
-#include "libSphere.h"
-#include "fuzzysphere.h"
-#include <stdio.h>
-#include <stdlib.h>
-
-// scattering from a uniform sphere w/ fuzzy surface
-// Modified from FuzzySpheres in libigor/libSphere.c without polydispersion: JHC
-double fuzzysphere_kernel(double dp[], double q){
-	double pi,x,xr;
-	double radius,sldSph,sldSolv,scale,bkg,delrho,fuzziness,f2,bes,vol,f;		//my local names
-
-	pi = 4.0*atan(1.0);
-	x= q;
-	scale = dp[0];
-	radius = dp[1];
-	fuzziness = dp[2];
-	sldSph = dp[3];
-	sldSolv = dp[4];
-	bkg = dp[5];
-	delrho=sldSph-sldSolv;
-
-	xr = x*radius;
-	//handle xr==0 separately
-	if(xr == 0.0){
-		bes = 1.0;
-	}else{
-		bes = 3.0*(sin(xr)-xr*cos(xr))/(xr*xr*xr);
-	}
-	vol = 4.0*pi/3.0*radius*radius*radius;
-	f = vol*bes*delrho;		// [=] A
-	f *= exp(-0.5*fuzziness*fuzziness*x*x);
-	// normalize to single particle volume, convert to 1/cm
-	f2 = f * f / vol * 1.0e8;		// [=] 1/cm
-
-	f2 *= scale;
-	f2 += bkg;
-
-    return(f2);	//scale, and add in the background
-}
-
-
-/**
- * Function to evaluate 1D scattering function
- * @param pars: parameters of the fuzzysphere
- * @param q: q-value
- * @return: function value
- */
-double fuzzysphere_analytical_1D(FuzzySphereParameters *pars, double q) {
-	double dp[6];
-
-	dp[0] = pars->scale;
-	dp[1] = pars->radius;
-	dp[2] = pars->fuzziness;
-	dp[3] = pars->sldSph;
-	dp[4] = pars->sldSolv;
-	dp[5] = pars->background;
-
-	return fuzzysphere_kernel(dp, q);
-}
-
-/**
- * Function to evaluate 2D scattering function
- * @param pars: parameters of the fuzzysphere
- * @param q: q-value
- * @return: function value
- */
-double fuzzysphere_analytical_2D(FuzzySphereParameters *pars, double q, double phi) {
-	return fuzzysphere_analytical_1D(pars,q);
-}
-
-double sphere_analytical_2DXY(FuzzySphereParameters *pars, double qx, double qy){
-	return fuzzysphere_analytical_1D(pars,sqrt(qx*qx+qy*qy));
-}
Index: sansmodels/src/c_extensions/fuzzysphere.h
===================================================================
--- sansmodels/src/c_extensions/fuzzysphere.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/fuzzysphere.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,59 +1,61 @@
 #if !defined(fuzzysphere_h)
 #define fuzzysphere_h
+#include "parameters.hh"
 
 /**
  * Structure definition for FuzzySphere parameters
  */
- //[PYTHONCLASS] = FuzzySphereModel
- //[DISP_PARAMS] = radius, fuzziness
- //[DESCRIPTION] =<text>
- //				scale: scale factor times volume fraction,
- //					or just volume fraction for absolute scale data
- //				radius: radius of the solid sphere
- //				fuzziness = the STD of the height of fuzzy interfacial
- //				 thickness (ie., so-called interfacial roughness)
- //				sldSph: the SLD of the sphere
- //				sldSolv: the SLD of the solvent
- //				background: incoherent background
- //			Note: By definition, this function works only when fuzziness << radius.
- //		</text>
- //[FIXED]=  radius.width; fuzziness.width
- //[ORIENTATION_PARAMS]= <text> </text>
+//[PYTHONCLASS] = FuzzySphereModel
+//[DISP_PARAMS] = radius, fuzziness
+//[DESCRIPTION] =<text>
+//				scale: scale factor times volume fraction,
+//					or just volume fraction for absolute scale data
+//				radius: radius of the solid sphere
+//				fuzziness = the STD of the height of fuzzy interfacial
+//				 thickness (ie., so-called interfacial roughness)
+//				sldSph: the SLD of the sphere
+//				sldSolv: the SLD of the solvent
+//				background: incoherent background
+//			Note: By definition, this function works only when fuzziness << radius.
+//		</text>
+//[FIXED]=  radius.width; fuzziness.width
+//[ORIENTATION_PARAMS]= <text> </text>
 
-typedef struct {
-    /// Scale factor
-    //  [DEFAULT]=scale= 0.01
-    double scale;
+class FuzzySphereModel{
+public:
+  // Model parameters
+  /// Radius of sphere [A]
+  //  [DEFAULT]=radius=60.0 [A]
+  Parameter radius;
+  /// Scale factor
+  //  [DEFAULT]=scale= 0.01
+  Parameter scale;
 
-    ///	Radius of sphere [A]
-    //  [DEFAULT]=radius=60.0 [A]
-    double radius;
 
-    ///	surface roughness [A]
-    //  [DEFAULT]=fuzziness= 10.0 [A]
-    double fuzziness;
+  /// surface roughness [A]
+  //  [DEFAULT]=fuzziness= 10.0 [A]
+  Parameter fuzziness;
 
-    ///	sldSph [1/A^(2)]
-    //  [DEFAULT]=sldSph= 1.0e-6 [1/A^(2)]
-    double sldSph;
+  /// sldSph [1/A^(2)]
+  //  [DEFAULT]=sldSph= 1.0e-6 [1/A^(2)]
+  Parameter sldSph;
 
-    ///	sldSolv [1/A^(2)]
-    //  [DEFAULT]=sldSolv= 3.0e-6 [1/A^(2)]
-    double sldSolv;
+  /// sldSolv [1/A^(2)]
+  //  [DEFAULT]=sldSolv= 3.0e-6 [1/A^(2)]
+  Parameter sldSolv;
 
-	/// Incoherent Background [1/cm]
-	//  [DEFAULT]=background=0.001 [1/cm]
-	double background;
-} FuzzySphereParameters;
+  /// Incoherent Background [1/cm]
+  //  [DEFAULT]=background=0.001 [1/cm]
+  Parameter background;
 
-/// kernel
-double fuzzysphere_kernel(double dp[], double q);
+  // Constructor
+  FuzzySphereModel();
 
-/// 1D scattering function
-double fuzzysphere_analytical_1D(FuzzySphereParameters *pars, double q);
+  // 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);
+};
 
-/// 2D scattering function
-double fuzzysphere_analytical_2D(FuzzySphereParameters *pars, double q, double phi);
-double fuzzysphere_analytical_2DXY(FuzzySphereParameters *pars, double qx, double qy);
-double fuzzysphere_analytical_2D_scaled(FuzzySphereParameters *pars, double q, double q_x, double q_y);
 #endif
Index: sansmodels/src/c_extensions/pearlnecklace.c
===================================================================
--- sansmodels/src/c_extensions/pearlnecklace.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,153 +1,0 @@
-/**
- * pearl necklace model
- */
-#include <math.h>
-#include "pearlnecklace.h"
-#include "libmultifunc/libfunc.h"
-#include <stdio.h>
-#include <stdlib.h>
-
-
-double pearl_necklace_kernel(double dp[], double q) {
-	// fit parameters
-	double scale = dp[0];
-	double radius = dp[1];
-	double edge_separation = dp[2];
-	double thick_string = dp[3];
-	double num_pearls = dp[4];
-	double sld_pearl = dp[5];
-	double sld_string = dp[6];
-	double sld_solv = dp[7];
-	double background = dp[8];
-
-	//relative slds
-	double contrast_pearl = sld_pearl - sld_solv;
-	double contrast_string = sld_string - sld_solv;
-
-	// number of string segments
-	double num_strings = num_pearls - 1.0;
-
-	//Pi
-	double pi = 4.0*atan(1.0);
-
-	// each volumes
-	double string_vol = edge_separation * pi * pow((thick_string / 2.0), 2);
-	double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3);
-
-	//total volume
-	double tot_vol;
-  //each masses
-  double m_r= contrast_string * string_vol;
-  double m_s= contrast_pearl * pearl_vol;
-  double psi;
-  double gamma;
-  double beta;
-  //form factors
-  double sss; //pearls
-  double srr; //strings
-  double srs; //cross
-  double A_s;
-  double srr_1;
-  double srr_2;
-  double srr_3;
-  double form_factor;
-
-  tot_vol = num_strings * string_vol;
-	tot_vol += num_pearls * pearl_vol;
-
-	//sine functions of a pearl
-	psi = sin(q*radius);
-	psi -= q * radius * cos(q * radius);
-	psi /= pow((q * radius), 3);
-	psi *= 3.0;
-
-	// Note take only 20 terms in Si series: 10 terms may be enough though.
-	gamma = Si(q* edge_separation);
-	gamma /= (q* edge_separation);
-	beta = Si(q * (edge_separation + radius));
-	beta -= Si(q * radius);
-	beta /= (q* edge_separation);
-
-	// center to center distance between the neighboring pearls
-	A_s = edge_separation + 2.0 * radius;
-
-	// form factor for num_pearls
-	sss = 1.0 - pow(sinc(q*A_s), num_pearls );
-	sss /= pow((1.0-sinc(q*A_s)), 2);
-	sss *= -sinc(q*A_s);
-	sss -= num_pearls/2.0;
-	sss += num_pearls/(1.0-sinc(q*A_s));
-	sss *= 2.0 * pow((m_s*psi), 2);
-
-	// form factor for num_strings (like thin rods)
-	srr_1 = -pow(sinc(q*edge_separation/2.0), 2);
-
-	srr_1 += 2.0 * gamma;
-	srr_1 *= num_strings;
-	srr_2 = 2.0/(1.0-sinc(q*A_s));
-	srr_2 *= num_strings;
-	srr_2 *= pow(beta, 2);
-	srr_3 = 1.0 - pow(sinc(q*A_s), num_strings);
-	srr_3 /= pow((1.0-sinc(q*A_s)), 2);
-	srr_3 *= pow(beta, 2);
-	srr_3 *= -2.0;
-
-	// total srr
-	srr = srr_1 + srr_2 + srr_3;
-	srr *= pow(m_r, 2);
-
-	// form factor for correlations
-	srs = 1.0;
-	srs -= pow(sinc(q*A_s), num_strings);
-	srs /= pow((1.0-sinc(q*A_s)), 2);
-	srs *= -sinc(q*A_s);
-	srs += (num_strings/(1.0-sinc(q*A_s)));
-	srs *= 4.0;
-	srs *= (m_r * m_s * beta * psi);
-
-	form_factor = sss + srr + srs;
-	form_factor /= (tot_vol * 1.0e-8); // norm by volume and A^-1 to cm^-1
-
-	// scale and background
-	form_factor *= scale;
-	form_factor += background;
-    return (form_factor);
-}
-
-/**
- * Function to evaluate pearlnecklace function
- * @param pars: parameters of pearlnecklace
- * @param q: q-value
- * @return: function value
- */
-
-double pearl_necklace_analytical_1D(PeralNecklaceParameters *pars, double q) {
-	double dp[9];
-
-	dp[0] = pars->scale;
-	dp[1] = pars->radius;
-	dp[2] = pars->edge_separation;
-	dp[3] = pars->thick_string;
-	dp[4] = pars->num_pearls;
-	dp[5] = pars->sld_pearl;
-	dp[6] = pars->sld_string;
-	dp[7] = pars->sld_solv;
-	dp[8] = pars->background;
-
-	return pearl_necklace_kernel(dp, q);
-}
-
-/**
- * Function to evaluate pearl_necklace function
- * @param pars: parameters of PeralNecklace
- * @param q: q-value
- * @return: function value
- */
-double pearl_necklace_analytical_2D(PeralNecklaceParameters *pars, double q, double phi) {
-	return pearl_necklace_analytical_1D(pars,q);
-}
-
-double pearl_necklace_analytical_2DXY(PeralNecklaceParameters *pars, double qx, double qy){
-	return pearl_necklace_analytical_1D(pars,sqrt(qx*qx+qy*qy));
-}
-
Index: sansmodels/src/c_extensions/pearlnecklace.h
===================================================================
--- sansmodels/src/c_extensions/pearlnecklace.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/pearlnecklace.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,66 +1,101 @@
 #if !defined(o_h)
 #define pearl_necklace_h
+#include "parameters.hh"
 
 /**
  * Structure definition for sphere parameters
  */
- //[PYTHONCLASS] = PearlNecklaceModel
- //[DISP_PARAMS] = radius, edge_separation
- //[DESCRIPTION] =<text>Calculate form factor for Pearl Necklace Model
+//[PYTHONCLASS] = PearlNecklaceModel
+//[DISP_PARAMS] = radius, edge_separation
+//[DESCRIPTION] =<text>Calculate form factor for Pearl Necklace Model
 //				[Macromol. Symp. 2004, 211, 25-42]
- //				Parameters:
- //				background:background
- //				scale: scale factor
- //				sld_pearl: the SLD of the pearl spheres
- //	            sld_string: the SLD of the strings
- //				sld_solv: the SLD of the solvent
- //				num_pearls: number of the pearls
- //				radius: the radius of a pearl
- //				edge_separation: the length of string segment; surface to surface
- //				thick_string: thickness (ie, diameter) of the string
- //		</text>
- //[FIXED]=  <text>radius.width; edge_separation.width</text>
- //[NON_FITTABLE_PARAMS]= <text></text>
- //[ORIENTATION_PARAMS]= <text> </text>
+//				Parameters:
+//				background:background
+//				scale: scale factor
+//				sld_pearl: the SLD of the pearl spheres
+//	            sld_string: the SLD of the strings
+//				sld_solv: the SLD of the solvent
+//				num_pearls: number of the pearls
+//				radius: the radius of a pearl
+//				edge_separation: the length of string segment; surface to surface
+//				thick_string: thickness (ie, diameter) of the string
+//		</text>
+//[FIXED]=  <text>radius.width; edge_separation.width</text>
+//[NON_FITTABLE_PARAMS]= <text></text>
+//[ORIENTATION_PARAMS]= <text> </text>
 
 typedef struct {
-    /// Scale factor
-    //  [DEFAULT]=scale= 1.0
-	double scale;
-    ///	radius [A]
-    //  [DEFAULT]=radius=80.0 [A]
-	double radius;
-	///	edge_separation
-	//  [DEFAULT]=edge_separation= 350 [A]
-	double edge_separation;
-	///	thick_string [A]
-	//  [DEFAULT]=thick_string= 2.5 [A]
-	double thick_string;
-	///	num_pearls
-	//  [DEFAULT]=num_pearls= 3
-	double num_pearls;
-	///	sld_pearl
-	//  [DEFAULT]=sld_pearl= 1.0e-06 [1/A^(2)]
-	double sld_pearl;
-	///	sld_string
-	//  [DEFAULT]=sld_string= 1.0e-06 [1/A^(2)]
-	double sld_string;
-	///	sld_solv
-	//  [DEFAULT]=sld_solv= 6.3e-06 [1/A^(2)]
-	double sld_solv;
-	/// Background
-	//  [DEFAULT]=background=0
-	double background;
+  /// Scale factor
+  //  [DEFAULT]=scale= 1.0
+  double scale;
+  ///	radius [A]
+  //  [DEFAULT]=radius=80.0 [A]
+  double radius;
+  ///	edge_separation
+  //  [DEFAULT]=edge_separation= 350 [A]
+  double edge_separation;
+  ///	thick_string [A]
+  //  [DEFAULT]=thick_string= 2.5 [A]
+  double thick_string;
+  ///	num_pearls
+  //  [DEFAULT]=num_pearls= 3
+  double num_pearls;
+  ///	sld_pearl
+  //  [DEFAULT]=sld_pearl= 1.0e-06 [1/A^(2)]
+  double sld_pearl;
+  ///	sld_string
+  //  [DEFAULT]=sld_string= 1.0e-06 [1/A^(2)]
+  double sld_string;
+  ///	sld_solv
+  //  [DEFAULT]=sld_solv= 6.3e-06 [1/A^(2)]
+  double sld_solv;
+  /// Background
+  //  [DEFAULT]=background=0
+  double background;
 
 }PeralNecklaceParameters;
 
-double pearl_necklace_kernel(double dq[], double q);
 
-/// 1D scattering function
-double pearl_necklace_analytical_1D(PeralNecklaceParameters *pars, double q);
 
-/// 2D scattering function
-double pearl_necklace_analytical_2D(PeralNecklaceParameters *pars, double q, double phi);
-double pearl_necklace_analytical_2DXY(PeralNecklaceParameters *pars, double qx, double qy);
+class PearlNecklaceModel{
+public:
+  // Model parameters
+  /// Scale factor
+  //  [DEFAULT]=scale= 1.0
+  Parameter scale;
+  /// radius [A]
+  //  [DEFAULT]=radius=80.0 [A]
+  Parameter radius;
+  /// edge_separation
+  //  [DEFAULT]=edge_separation= 350 [A]
+  Parameter edge_separation;
+  /// thick_string [A]
+  //  [DEFAULT]=thick_string= 2.5 [A]
+  Parameter thick_string;
+  /// num_pearls
+  //  [DEFAULT]=num_pearls= 3
+  Parameter num_pearls;
+  /// sld_pearl
+  //  [DEFAULT]=sld_pearl= 1.0e-06 [1/A^(2)]
+  Parameter sld_pearl;
+  /// sld_string
+  //  [DEFAULT]=sld_string= 1.0e-06 [1/A^(2)]
+  Parameter sld_string;
+  /// sld_solv
+  //  [DEFAULT]=sld_solv= 6.3e-06 [1/A^(2)]
+  Parameter sld_solv;
+  /// Background
+  //  [DEFAULT]=background=0
+  Parameter background;
+
+  // Constructor
+  PearlNecklaceModel();
+
+  // 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_extensions/sld_cal.c
===================================================================
--- sansmodels/src/c_extensions/sld_cal.c	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ 	(revision )
@@ -1,47 +1,0 @@
-/**
- * SLD Calculator
- * @author: Jae Hie Cho / UTK
- */
-
-#include "sld_cal.h"
-#include "libmultifunc/librefl.h"
-#include <math.h>
-
-
-/**
- * Function to calculate sld
- * @param pars: parameters
- * @param x: independent param-value
- * @return: sld value
- */
-double sld_cal_analytical_1D(SLDCalParameters *pars, double x) {
-	double fun, nsl, n_s, fun_coef, sld_l, sld_r, sld_out;
-	int fun_type;
-
-	fun = pars->fun_type;
-	nsl = pars->npts_inter;
-	n_s = pars->shell_num;
-	fun_coef = pars->nu_inter;
-	sld_l = pars-> sld_left;
-	sld_r = pars-> sld_right;
-
-	fun_type = floor(fun);
-
-	sld_out = intersldfunc(fun_type, nsl, n_s, fun_coef, sld_l, sld_r);
-
-	return sld_out;
-}
-
-/**
- * Function to evaluate 2D scattering function
- * @param pars: parameters
- * @param q: q-value
- * @return: function value
- */
-double sld_cal_analytical_2D(SLDCalParameters *pars, double x, double phi) {
-	return sld_cal_analytical_1D(pars,x);
-}
-
-double sld_cal_analytical_2DXY(SLDCalParameters *pars, double xx, double xy){
-	return sld_cal_analytical_1D(pars,sqrt(xx*xx+xy*xy));
-}
Index: sansmodels/src/c_extensions/sld_cal.h
===================================================================
--- sansmodels/src/c_extensions/sld_cal.h	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_extensions/sld_cal.h	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,4 +1,5 @@
 #if !defined(sld_cal_h)
 #define sld_cal_h
+#include "parameters.hh"
 
 /**
@@ -12,39 +13,43 @@
 
  **/
-typedef struct {
-    ///	fun_type
-    //  [DEFAULT]=fun_type=0
-    double fun_type;
 
-    /// npts_inter
-    //  [DEFAULT]=npts_inter= 21
-    double npts_inter;
+class SLDCalFunc{
+public:
+  // Model parameters
+  /// fun_type
+  //  [DEFAULT]=fun_type=0
+  Parameter fun_type;
 
-    /// shell_num
-    //  [DEFAULT]=shell_num= 0
-    double shell_num;
+  /// npts_inter
+  //  [DEFAULT]=npts_inter= 21
+  Parameter npts_inter;
 
-    /// nu_inter
-    //  [DEFAULT]=nu_inter= 2.5
-    double nu_inter;
+  /// shell_num
+  //  [DEFAULT]=shell_num= 0
+  Parameter shell_num;
 
-    /// sld_left [1/A^(2)]
-    //  [DEFAULT]=sld_left= 0 [1/A^(2)]
-    double sld_left;
+  /// nu_inter
+  //  [DEFAULT]=nu_inter= 2.5
+  Parameter nu_inter;
 
-    /// sld_right [1/A^(2)]
-    //  [DEFAULT]=sld_right= 0 [1/A^(2)]
-    double sld_right;
+  /// sld_left [1/A^(2)]
+  //  [DEFAULT]=sld_left= 0 [1/A^(2)]
+  Parameter sld_left;
 
-} SLDCalParameters;
+  /// sld_right [1/A^(2)]
+  //  [DEFAULT]=sld_right= 0 [1/A^(2)]
+  Parameter sld_right;
+
+  // Constructor
+  SLDCalFunc();
+
+  // Operators to get SLD
+  double operator()(double q);
+  double operator()(double qx, double qy);
+  double calculate_ER();
+  double evaluate_rphi(double q, double phi);
+};
 
 
 
-/// 1D function
-double sld_cal_analytical_1D(SLDCalParameters *pars, double q);
-
-/// 2D function
-double sld_cal_analytical_2D(SLDCalParameters *pars, double q, double phi);
-double sld_cal_analytical_2DXY(SLDCalParameters *pars, double qx, double qy);
-
 #endif
Index: sansmodels/src/c_models/DiamCyl.cpp
===================================================================
--- sansmodels/src/c_models/DiamCyl.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/DiamCyl.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -18,16 +18,14 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "DiamCyl.h"
 
 extern "C" {
 	#include "libStructureFactor.h"
-	#include "DiamCyl.h"
 }
 
@@ -112,37 +110,2 @@
   return 0.0;
 }
-// Testing code
-
-/**
-int main(void)
-{
-	DiamCylFunc c = DiamCylFunc();
-
-	printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	c.radius.dispersion = new GaussianDispersion();
-	c.radius.dispersion->npts = 100;
-	c.radius.dispersion->width = 20;
-
-	//c.length.dispersion = GaussianDispersion();
-	//c.length.dispersion.npts = 20;
-	//c.length.dispersion.width = 65;
-
-	//printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	//printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	//printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
-	//printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854));
-
-
-
-	double i_avg = c(0.01, 0.01);
-	double i_1d = c(sqrt(0.0002));
-
-	printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg);
-	printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d);
-	printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg);
-
-
-	return 0;
-}
- **/
Index: sansmodels/src/c_models/DiamEllip.cpp
===================================================================
--- sansmodels/src/c_models/DiamEllip.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/DiamEllip.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -21,12 +21,11 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "DiamEllip.h"
 
 extern "C" {
 	#include "libStructureFactor.h"
-	#include "DiamEllip.h"
 }
 
@@ -109,36 +108,2 @@
   return 0.0;
 }
-// Testing code
-/*
-int main(void)
-{
-	SquareWellModel c = SquareWellModel();
-
-	printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	c.radius.dispersion = new GaussianDispersion();
-	c.radius.dispersion->npts = 100;
-	c.radius.dispersion->width = 5;
-
-	//c.length.dispersion = GaussianDispersion();
-	//c.length.dispersion.npts = 20;
-	//c.length.dispersion.width = 65;
-
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
-	printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854));
-
-
-
-	double i_avg = c(0.01, 0.01);
-	double i_1d = c(sqrt(0.0002));
-
-	printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg);
-	printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d);
-	printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg);
-
-
-	return 0;
-}
-*/
Index: sansmodels/src/c_models/Hardsphere.cpp
===================================================================
--- sansmodels/src/c_models/Hardsphere.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/Hardsphere.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -22,19 +22,18 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "Hardsphere.h"
 
 extern "C" {
-	#include "libStructureFactor.h"
-	#include "Hardsphere.h"
+#include "libStructureFactor.h"
 }
 
 HardsphereStructure :: HardsphereStructure() {
-	effect_radius      = Parameter(50.0, true);
-	effect_radius.set_min(0.0);
-	volfraction = Parameter(0.20, true);
-	volfraction.set_min(0.0);
+  effect_radius      = Parameter(50.0, true);
+  effect_radius.set_min(0.0);
+  volfraction = Parameter(0.20, true);
+  volfraction.set_min(0.0);
 }
 
@@ -46,28 +45,28 @@
  */
 double HardsphereStructure :: operator()(double q) {
-	double dp[2];
+  double dp[2];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = effect_radius();
-	dp[1] = volfraction();
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = effect_radius();
+  dp[1] = volfraction();
 
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	effect_radius.get_weights(weights_rad);
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  effect_radius.get_weights(weights_rad);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp[0] = weights_rad[i].value;
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[0] = weights_rad[i].value;
 
-		sum += weights_rad[i].weight
-				* HardSphereStruct(dp, q);
-		norm += weights_rad[i].weight;
-	}
-	return sum/norm ;
+    sum += weights_rad[i].weight
+        * HardSphereStruct(dp, q);
+    norm += weights_rad[i].weight;
+  }
+  return sum/norm ;
 }
 
@@ -79,6 +78,6 @@
  */
 double HardsphereStructure :: 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);
 }
 
@@ -91,7 +90,7 @@
  */
 double HardsphereStructure :: 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);
 }
 /**
@@ -100,5 +99,5 @@
  */
 double HardsphereStructure :: calculate_ER() {
-//NOT implemented yet!!!
+  //NOT implemented yet!!!
   return 0.0;
 }
Index: sansmodels/src/c_models/HayterMSA.cpp
===================================================================
--- sansmodels/src/c_models/HayterMSA.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/HayterMSA.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -22,24 +22,23 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "HayterMSA.h"
 
 extern "C" {
-	#include "libStructureFactor.h"
-	#include "HayterMSA.h"
+#include "libStructureFactor.h"
 }
 
 HayterMSAStructure :: HayterMSAStructure() {
-	effect_radius      = Parameter(20.75, true);
-	effect_radius.set_min(0.0);
-	charge      = Parameter(19.0, true);
-	volfraction = Parameter(0.0192, true);
-	volfraction.set_min(0.0);
-	temperature = Parameter(318.16, true);
-	temperature.set_min(0.0);
-	saltconc   = Parameter(0.0);
-	dielectconst  = Parameter(71.08);
+  effect_radius      = Parameter(20.75, true);
+  effect_radius.set_min(0.0);
+  charge      = Parameter(19.0, true);
+  volfraction = Parameter(0.0192, true);
+  volfraction.set_min(0.0);
+  temperature = Parameter(318.16, true);
+  temperature.set_min(0.0);
+  saltconc   = Parameter(0.0);
+  dielectconst  = Parameter(71.08);
 }
 
@@ -51,32 +50,32 @@
  */
 double HayterMSAStructure :: operator()(double q) {
-	double dp[6];
+  double dp[6];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = 2.0*effect_radius();
-	dp[1] = fabs(charge());
-	dp[2] = volfraction();
-	dp[3] = temperature();
-	dp[4] = saltconc();
-	dp[5] = dielectconst();
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = 2.0*effect_radius();
+  dp[1] = fabs(charge());
+  dp[2] = volfraction();
+  dp[3] = temperature();
+  dp[4] = saltconc();
+  dp[5] = dielectconst();
 
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	effect_radius.get_weights(weights_rad);
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  effect_radius.get_weights(weights_rad);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp[0] = 2.0*weights_rad[i].value;
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[0] = 2.0*weights_rad[i].value;
 
-		sum += weights_rad[i].weight
-				* HayterPenfoldMSA(dp, q);
-		norm += weights_rad[i].weight;
-	}
-	return sum/norm ;
+    sum += weights_rad[i].weight
+        * HayterPenfoldMSA(dp, q);
+    norm += weights_rad[i].weight;
+  }
+  return sum/norm ;
 }
 
@@ -88,6 +87,6 @@
  */
 double HayterMSAStructure :: 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);
 }
 
@@ -100,7 +99,7 @@
  */
 double HayterMSAStructure :: 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);
 }
 /**
@@ -109,40 +108,5 @@
  */
 double HayterMSAStructure :: calculate_ER() {
-//NOT implemented yet!!!
+  //NOT implemented yet!!!
   return 0.0;
 }
-
-// Testing code
-/*
-int main(void)
-{
-	SquareWellModel c = SquareWellModel();
-
-	printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	c.radius.dispersion = new GaussianDispersion();
-	c.radius.dispersion->npts = 100;
-	c.radius.dispersion->width = 5;
-
-	//c.length.dispersion = GaussianDispersion();
-	//c.length.dispersion.npts = 20;
-	//c.length.dispersion.width = 65;
-
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
-	printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854));
-
-
-
-	double i_avg = c(0.01, 0.01);
-	double i_1d = c(sqrt(0.0002));
-
-	printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg);
-	printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d);
-	printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg);
-
-
-	return 0;
-}
-*/
Index: sansmodels/src/c_models/SquareWell.cpp
===================================================================
--- sansmodels/src/c_models/SquareWell.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/SquareWell.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -18,25 +18,23 @@
  *   sansmodels/src/libigor
  *
- *	TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "SquareWell.h"
 
 extern "C" {
-	#include "libStructureFactor.h"
-	#include "SquareWell.h"
+#include "libStructureFactor.h"
 }
 
 SquareWellStructure :: SquareWellStructure() {
-	effect_radius      = Parameter(50.0, true);
-	effect_radius.set_min(0.0);
-	volfraction = Parameter(0.04, true);
-	volfraction.set_min(0.0);
-	welldepth   = Parameter(1.50);
-	wellwidth  = Parameter(1.20);
+  effect_radius      = Parameter(50.0, true);
+  effect_radius.set_min(0.0);
+  volfraction = Parameter(0.04, true);
+  volfraction.set_min(0.0);
+  welldepth   = Parameter(1.50);
+  wellwidth  = Parameter(1.20);
 }
 
@@ -48,30 +46,30 @@
  */
 double SquareWellStructure :: operator()(double q) {
-	double dp[4];
+  double dp[4];
 
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = effect_radius();
-	dp[1] = volfraction();
-	dp[2] = welldepth();
-	dp[3] = wellwidth();
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = effect_radius();
+  dp[1] = volfraction();
+  dp[2] = welldepth();
+  dp[3] = wellwidth();
 
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_rad;
-	effect_radius.get_weights(weights_rad);
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  effect_radius.get_weights(weights_rad);
 
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-	// Loop over radius weight points
-	for(size_t i=0; i<weights_rad.size(); i++) {
-		dp[0] = weights_rad[i].value;
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[0] = weights_rad[i].value;
 
-		sum += weights_rad[i].weight
-				* SquareWellStruct(dp, q);
-		norm += weights_rad[i].weight;
-	}
-	return sum/norm ;
+    sum += weights_rad[i].weight
+        * SquareWellStruct(dp, q);
+    norm += weights_rad[i].weight;
+  }
+  return sum/norm ;
 }
 
@@ -83,6 +81,6 @@
  */
 double SquareWellStructure :: 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);
 }
 
@@ -95,7 +93,7 @@
  */
 double SquareWellStructure :: 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);
 }
 /**
@@ -104,39 +102,5 @@
  */
 double SquareWellStructure :: calculate_ER() {
-//NOT implemented yet!!!
+  //NOT implemented yet!!!
   return 0.0;
 }
-// Testing code
-/*
-int main(void)
-{
-	SquareWellModel c = SquareWellModel();
-
-	printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	c.radius.dispersion = new GaussianDispersion();
-	c.radius.dispersion->npts = 100;
-	c.radius.dispersion->width = 5;
-
-	//c.length.dispersion = GaussianDispersion();
-	//c.length.dispersion.npts = 20;
-	//c.length.dispersion.width = 65;
-
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	printf("I(Q=%g) = %g\n", 0.001, c(0.001));
-	printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
-	printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854));
-
-
-
-	double i_avg = c(0.01, 0.01);
-	double i_1d = c(sqrt(0.0002));
-
-	printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg);
-	printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d);
-	printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg);
-
-
-	return 0;
-}
-*/
Index: sansmodels/src/c_models/StickyHS.cpp
===================================================================
--- sansmodels/src/c_models/StickyHS.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/StickyHS.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -22,12 +22,11 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "StickyHS.h"
 
 extern "C" {
 	#include "libStructureFactor.h"
-	#include "StickyHS.h"
 }
 
Index: sansmodels/src/c_models/bcc.cpp
===================================================================
--- sansmodels/src/c_models/bcc.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/bcc.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -21,12 +21,134 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "bcc.h"
 
 extern "C" {
 	#include "libSphere.h"
-	#include "bcc.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;
+
+} BCParameters;
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the BCCCrystalModel
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double bc_analytical_2D_scaled(BCParameters *pars, double q, double q_x, double q_y) {
+  double b3_x, b3_y, b3_z, b1_x, b1_y;
+  double q_z;
+  double alpha, cos_val_b3, cos_val_b2, cos_val_b1;
+  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, Fkq, Fkq_2;
+
+  //convert angle degree to radian
+  double pi = 4.0*atan(1.0);
+  double theta = pars->theta * pi/180.0;
+  double phi = pars->phi * pi/180.0;
+  double psi = pars->psi * pi/180.0;
+
+  double dp[5];
+  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);
+
+  //the occupied volume of the lattice
+  latticeScale = 2.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa/sqrt(3.0/4.0),3.0);
+  // q vector
+  q_z = 0.0; // for SANS; assuming qz is negligible
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
+    // b3 axis orientation
+    b3_x = sin(theta) * cos(phi);//negative sign here???
+    b3_y = sin(theta) * sin(phi);
+    b3_z = cos(theta);
+    cos_val_b3 =  b3_x*q_x + b3_y*q_y + b3_z*q_z;
+
+    alpha = acos(cos_val_b3);
+    // b1 axis orientation
+    b1_x = sin(psi);
+    b1_y = cos(psi);
+  cos_val_b1 = (b1_x*q_x + b1_y*q_y);
+    // b2 axis orientation
+  cos_val_b2 = sin(acos(cos_val_b1));
+  // alpha corrections
+  cos_val_b2 *= sin(alpha);
+  cos_val_b1 *= sin(alpha);
+
+    // Compute the angle btw vector q and the a3 axis
+    a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1-cos_val_b3);
+
+    // a1 axis
+    a1_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b2-cos_val_b1);
+
+    // a2 axis
+    a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1-cos_val_b2);
+
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      printf("bcc_ana_2D: Unexpected error: cos(alpha)>1\n");
+      return 0;
+    }
+    // Get Fkq and Fkq_2
+    Fkq = exp(-0.5*pow(Da/aa,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+    Fkq_2 = Fkq*Fkq;
+    // Call Zq=Z1*Z2*Z3
+    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_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 BCC_ParaCrystal
+ * @param q: q-value
+ * @return: function value
+ */
+static double bc_analytical_2DXY(BCParameters *pars, double qx, double qy){
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return bc_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
Index: sansmodels/src/c_models/fcc.cpp
===================================================================
--- sansmodels/src/c_models/fcc.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/fcc.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -21,26 +21,144 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "fcc.h"
 
 extern "C" {
-	#include "libSphere.h"
-	#include "fcc.h"
+#include "libSphere.h"
+}
+
+// Convenience parameter structure
+typedef struct {
+  double scale;
+  double dnn;
+  double d_factor;
+  double radius;
+  double sldSph;
+  double sldSolv;
+  double background;
+  double theta;
+  double phi;
+  double psi;
+} FCParameters;
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the FCCCrystalModel
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double fc_analytical_2D_scaled(FCParameters *pars, double q, double q_x, double q_y) {
+  double b3_x, b3_y, b3_z, b1_x, b1_y;
+  double q_z;
+  double alpha, cos_val_b3, cos_val_b2, cos_val_b1;
+  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, Fkq, Fkq_2;
+
+  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);
+  //contrast = pars->sldSph - pars->sldSolv;
+
+  latticeScale = 4.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa*sqrt(2.0),3.0);
+  // q vector
+  q_z = 0.0; // for SANS; assuming qz is negligible
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate in PRB 36(46), 3(6), 1754(3854)
+  // b3 axis orientation
+  b3_x = sin(theta) * cos(phi);//negative sign here???
+  b3_y = sin(theta) * sin(phi);
+  b3_z = cos(theta);
+  cos_val_b3 =  b3_x*q_x + b3_y*q_y + b3_z*q_z;
+  alpha = acos(cos_val_b3);
+  // b1 axis orientation
+  b1_x = sin(psi);
+  b1_y = cos(psi);
+  cos_val_b1 = (b1_x*q_x + b1_y*q_y);
+  // b2 axis orientation
+  cos_val_b2 = sin(acos(cos_val_b1));
+  // alpha correction
+  cos_val_b2 *= sin(alpha);
+  cos_val_b1 *= sin(alpha);
+
+  // Compute the angle btw vector q and the a3 axis
+  a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1);
+
+  // a1 axis
+  a1_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b3);
+
+  // a2 axis
+  a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1);
+
+  // The following test should always pass
+  if (fabs(cos_val_b3)>1.0) {
+    printf("fcc_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+  // Get Fkq and Fkq_2
+  Fkq = exp(-0.5*pow(Da/aa,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+  Fkq_2 = Fkq*Fkq;
+  // Call Zq=Z1*Z2*Z3
+  Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+  Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+  Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_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 FCC_ParaCrystal
+ * @param q: q-value
+ * @return: function value
+ */
+static double fc_analytical_2DXY(FCParameters *pars, double qx, double qy){
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return fc_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 FCCrystalModel :: FCCrystalModel() {
-	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 +170,49 @@
  */
 double FCCrystalModel :: 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 = FCC_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 = FCC_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,84 +224,84 @@
  */
 double FCCrystalModel :: operator()(double qx, double qy) {
-	FCParameters 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
-										* fc_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!!!
-					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();
+  FCParameters 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
+              * fc_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!!!
+          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();
 }
 
@@ -196,5 +314,5 @@
  */
 double FCCrystalModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -204,5 +322,5 @@
  */
 double FCCrystalModel :: calculate_ER() {
-	//NOT implemented yet!!!
-	return 0.0;
-}
+  //NOT implemented yet!!!
+  return 0.0;
+}
Index: sansmodels/src/c_models/fuzzysphere.cpp
===================================================================
--- sansmodels/src/c_models/fuzzysphere.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/fuzzysphere.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -15,11 +15,47 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
+#include "fuzzysphere.h"
 
 extern "C" {
-	#include "fuzzysphere.h"
+#include "libSphere.h"
+}
+
+// scattering from a uniform sphere w/ fuzzy surface
+// Modified from FuzzySpheres in libigor/libSphere.c without polydispersion: JHC
+static double fuzzysphere_kernel(double dp[], double q){
+  double pi,x,xr;
+  double radius,sldSph,sldSolv,scale,bkg,delrho,fuzziness,f2,bes,vol,f;   //my local names
+
+  pi = 4.0*atan(1.0);
+  x= q;
+  scale = dp[0];
+  radius = dp[1];
+  fuzziness = dp[2];
+  sldSph = dp[3];
+  sldSolv = dp[4];
+  bkg = dp[5];
+  delrho=sldSph-sldSolv;
+
+  xr = x*radius;
+  //handle xr==0 separately
+  if(xr == 0.0){
+    bes = 1.0;
+  }else{
+    bes = 3.0*(sin(xr)-xr*cos(xr))/(xr*xr*xr);
+  }
+  vol = 4.0*pi/3.0*radius*radius*radius;
+  f = vol*bes*delrho;   // [=] A
+  f *= exp(-0.5*fuzziness*fuzziness*x*x);
+  // normalize to single particle volume, convert to 1/cm
+  f2 = f * f / vol * 1.0e8;   // [=] 1/cm
+
+  f2 *= scale;
+  f2 += bkg;
+
+    return(f2); //scale, and add in the background
 }
 
Index: sansmodels/src/c_models/models.hh
===================================================================
--- sansmodels/src/c_models/models.hh	(revision 6e10cffb542a4045bc3c41ac60c4abf3ed220bb2)
+++ sansmodels/src/c_models/models.hh	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -27,223 +27,4 @@
 using namespace std;
 
-
-
-class PearlNecklaceModel{
-public:
-	// Model parameters
-	Parameter scale;
-	Parameter radius;
-	Parameter edge_separation;
-	Parameter thick_string;
-	Parameter num_pearls;
-	Parameter sld_pearl;
-	Parameter sld_string;
-	Parameter sld_solv;
-	Parameter background;
-
-	// Constructor
-	PearlNecklaceModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class FCCrystalModel{
-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
-	FCCrystalModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-
-class BCCrystalModel{
-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
-	BCCrystalModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-
-class FuzzySphereModel{
-public:
-	// Model parameters
-	Parameter radius;
-	Parameter scale;
-	Parameter fuzziness;
-	Parameter sldSph;
-	Parameter sldSolv;
-	Parameter background;
-
-	// Constructor
-	FuzzySphereModel();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class HardsphereStructure{
-public:
-	// Model parameters
-	Parameter effect_radius;
-	Parameter volfraction;
-
-	// Constructor
-	HardsphereStructure();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class StickyHSStructure{
-public:
-	// Model parameters
-	Parameter effect_radius;
-	Parameter volfraction;
-	Parameter perturb;
-	Parameter stickiness;
-
-	// Constructor
-	StickyHSStructure();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class SquareWellStructure{
-public:
-	// Model parameters
-	Parameter effect_radius;
-	Parameter volfraction;
-	Parameter welldepth;
-	Parameter wellwidth;
-
-	// Constructor
-	SquareWellStructure();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class HayterMSAStructure{
-public:
-	// Model parameters
-	Parameter effect_radius;
-	Parameter charge;
-	Parameter volfraction;
-	Parameter temperature;
-	Parameter saltconc;
-	Parameter dielectconst;
-
-	// Constructor
-	HayterMSAStructure();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class DiamEllipFunc{
-public:
-	// Model parameters
-	Parameter radius_a;
-	Parameter radius_b;
-
-	// Constructor
-	DiamEllipFunc();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-class DiamCylFunc{
-public:
-	// Model parameters
-	Parameter radius;
-	Parameter length;
-
-	// Constructor
-	DiamCylFunc();
-
-	// Operators to get I(Q)
-	double operator()(double q);
-	double operator()(double qx, double qy);
-	double calculate_ER();
-	double evaluate_rphi(double q, double phi);
-};
-
-
-class SLDCalFunc{
-public:
-	// Model parameters
-	Parameter fun_type;
-	Parameter npts_inter;
-	Parameter shell_num;
-	Parameter nu_inter;
-	Parameter sld_left;
-	Parameter sld_right;
-
-	// Constructor
-	SLDCalFunc();
-
-	// Operators to get SLD
-	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/pearlnecklace.cpp
===================================================================
--- sansmodels/src/c_models/pearlnecklace.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/pearlnecklace.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -1,27 +1,133 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "pearlnecklace.h"
 
 extern "C" {
-	#include "pearlnecklace.h"
+#include "libmultifunc/libfunc.h"
+}
+
+static double pearl_necklace_kernel(double dp[], double q) {
+  // fit parameters
+  double scale = dp[0];
+  double radius = dp[1];
+  double edge_separation = dp[2];
+  double thick_string = dp[3];
+  double num_pearls = dp[4];
+  double sld_pearl = dp[5];
+  double sld_string = dp[6];
+  double sld_solv = dp[7];
+  double background = dp[8];
+
+  //relative slds
+  double contrast_pearl = sld_pearl - sld_solv;
+  double contrast_string = sld_string - sld_solv;
+
+  // number of string segments
+  double num_strings = num_pearls - 1.0;
+
+  //Pi
+  double pi = 4.0*atan(1.0);
+
+  // each volumes
+  double string_vol = edge_separation * pi * pow((thick_string / 2.0), 2);
+  double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3);
+
+  //total volume
+  double tot_vol;
+  //each masses
+  double m_r= contrast_string * string_vol;
+  double m_s= contrast_pearl * pearl_vol;
+  double psi;
+  double gamma;
+  double beta;
+  //form factors
+  double sss; //pearls
+  double srr; //strings
+  double srs; //cross
+  double A_s;
+  double srr_1;
+  double srr_2;
+  double srr_3;
+  double form_factor;
+
+  tot_vol = num_strings * string_vol;
+  tot_vol += num_pearls * pearl_vol;
+
+  //sine functions of a pearl
+  psi = sin(q*radius);
+  psi -= q * radius * cos(q * radius);
+  psi /= pow((q * radius), 3);
+  psi *= 3.0;
+
+  // Note take only 20 terms in Si series: 10 terms may be enough though.
+  gamma = Si(q* edge_separation);
+  gamma /= (q* edge_separation);
+  beta = Si(q * (edge_separation + radius));
+  beta -= Si(q * radius);
+  beta /= (q* edge_separation);
+
+  // center to center distance between the neighboring pearls
+  A_s = edge_separation + 2.0 * radius;
+
+  // form factor for num_pearls
+  sss = 1.0 - pow(sinc(q*A_s), num_pearls );
+  sss /= pow((1.0-sinc(q*A_s)), 2);
+  sss *= -sinc(q*A_s);
+  sss -= num_pearls/2.0;
+  sss += num_pearls/(1.0-sinc(q*A_s));
+  sss *= 2.0 * pow((m_s*psi), 2);
+
+  // form factor for num_strings (like thin rods)
+  srr_1 = -pow(sinc(q*edge_separation/2.0), 2);
+
+  srr_1 += 2.0 * gamma;
+  srr_1 *= num_strings;
+  srr_2 = 2.0/(1.0-sinc(q*A_s));
+  srr_2 *= num_strings;
+  srr_2 *= pow(beta, 2);
+  srr_3 = 1.0 - pow(sinc(q*A_s), num_strings);
+  srr_3 /= pow((1.0-sinc(q*A_s)), 2);
+  srr_3 *= pow(beta, 2);
+  srr_3 *= -2.0;
+
+  // total srr
+  srr = srr_1 + srr_2 + srr_3;
+  srr *= pow(m_r, 2);
+
+  // form factor for correlations
+  srs = 1.0;
+  srs -= pow(sinc(q*A_s), num_strings);
+  srs /= pow((1.0-sinc(q*A_s)), 2);
+  srs *= -sinc(q*A_s);
+  srs += (num_strings/(1.0-sinc(q*A_s)));
+  srs *= 4.0;
+  srs *= (m_r * m_s * beta * psi);
+
+  form_factor = sss + srr + srs;
+  form_factor /= (tot_vol * 1.0e-8); // norm by volume and A^-1 to cm^-1
+
+  // scale and background
+  form_factor *= scale;
+  form_factor += background;
+  return (form_factor);
 }
 
 PearlNecklaceModel :: PearlNecklaceModel() {
-	scale = Parameter(1.0);
-	radius = Parameter(80.0, true);
-	radius.set_min(0.0);
-	edge_separation = Parameter(350.0, true);
-	edge_separation.set_min(0.0);
-	thick_string = Parameter(2.5, true);
-	thick_string.set_min(0.0);
-	num_pearls = Parameter(3);
-	num_pearls.set_min(0.0);
-	sld_pearl = Parameter(1.0e-06);
-	sld_string = Parameter(1.0e-06);
-	sld_solv = Parameter(6.3e-06);
-    background = Parameter(0.0);
+  scale = Parameter(1.0);
+  radius = Parameter(80.0, true);
+  radius.set_min(0.0);
+  edge_separation = Parameter(350.0, true);
+  edge_separation.set_min(0.0);
+  thick_string = Parameter(2.5, true);
+  thick_string.set_min(0.0);
+  num_pearls = Parameter(3);
+  num_pearls.set_min(0.0);
+  sld_pearl = Parameter(1.0e-06);
+  sld_string = Parameter(1.0e-06);
+  sld_solv = Parameter(6.3e-06);
+  background = Parameter(0.0);
 
 }
@@ -33,55 +139,55 @@
  */
 double PearlNecklaceModel :: operator()(double q) {
-	double dp[9];
-	// Fill parameter array for IGOR library
-	// Add the background after averaging
-	dp[0] = scale();
-	dp[1] = radius();
-	dp[2] = edge_separation();
-	dp[3] = thick_string();
-	dp[4] = num_pearls();
-	dp[5] = sld_pearl();
-	dp[6] = sld_string();
-	dp[7] = sld_solv();
-	dp[8] = 0.0;
-	double pi = 4.0*atan(1.0);
-	// No polydispersion supported in this model.
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_radius;
-	radius.get_weights(weights_radius);
-	vector<WeightPoint> weights_edge_separation;
-	edge_separation.get_weights(weights_edge_separation);
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double vol = 0.0;
-	double string_vol = 0.0;
-	double pearl_vol = 0.0;
-	double tot_vol = 0.0;
-	// Loop over core weight points
-	for(size_t i=0; i<weights_radius.size(); i++) {
-		dp[1] = weights_radius[i].value;
-		// Loop over thick_inter0 weight points
-		for(size_t j=0; j<weights_edge_separation.size(); j++) {
-			dp[2] = weights_edge_separation[j].value;
-			pearl_vol = 4.0 /3.0 * pi * pow(dp[1], 3);
-			string_vol =dp[2] * pi * pow((dp[3] / 2.0), 2);
-			tot_vol = (dp[4] - 1.0) * string_vol;
-			tot_vol += dp[4] * pearl_vol;
-			//Un-normalize Sphere by volume
-			sum += weights_radius[i].weight * weights_edge_separation[j].weight
-				* pearl_necklace_kernel(dp,q) * tot_vol;
-			//Find average volume
-			vol += weights_radius[i].weight * weights_edge_separation[j].weight
-				* tot_vol;
-			norm += weights_radius[i].weight * weights_edge_separation[j].weight;
-		}
-	}
-
-	if (vol != 0.0 && norm != 0.0) {
-		//Re-normalize by avg volume
-		sum = sum/(vol/norm);}
-
-	return sum/norm + background();
+  double dp[9];
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = radius();
+  dp[2] = edge_separation();
+  dp[3] = thick_string();
+  dp[4] = num_pearls();
+  dp[5] = sld_pearl();
+  dp[6] = sld_string();
+  dp[7] = sld_solv();
+  dp[8] = 0.0;
+  double pi = 4.0*atan(1.0);
+  // No polydispersion supported in this model.
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+  vector<WeightPoint> weights_edge_separation;
+  edge_separation.get_weights(weights_edge_separation);
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+  double string_vol = 0.0;
+  double pearl_vol = 0.0;
+  double tot_vol = 0.0;
+  // Loop over core weight points
+  for(size_t i=0; i<weights_radius.size(); i++) {
+    dp[1] = weights_radius[i].value;
+    // Loop over thick_inter0 weight points
+    for(size_t j=0; j<weights_edge_separation.size(); j++) {
+      dp[2] = weights_edge_separation[j].value;
+      pearl_vol = 4.0 /3.0 * pi * pow(dp[1], 3);
+      string_vol =dp[2] * pi * pow((dp[3] / 2.0), 2);
+      tot_vol = (dp[4] - 1.0) * string_vol;
+      tot_vol += dp[4] * pearl_vol;
+      //Un-normalize Sphere by volume
+      sum += weights_radius[i].weight * weights_edge_separation[j].weight
+          * pearl_necklace_kernel(dp,q) * tot_vol;
+      //Find average volume
+      vol += weights_radius[i].weight * weights_edge_separation[j].weight
+          * tot_vol;
+      norm += weights_radius[i].weight * weights_edge_separation[j].weight;
+    }
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -93,6 +199,6 @@
  */
 double PearlNecklaceModel :: 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);
 }
 
@@ -105,5 +211,5 @@
  */
 double PearlNecklaceModel :: evaluate_rphi(double q, double phi) {
-	return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -118,60 +224,60 @@
 // sld of solvent.
 double PearlNecklaceModel :: calculate_ER() {
-	PeralNecklaceParameters dp;
-
-	dp.scale = scale();
-	dp.radius = radius();
-	dp.edge_separation = edge_separation();
-	dp.thick_string = thick_string();
-	dp.num_pearls = num_pearls();
-
-	double rad_out = 0.0;
-	// Perform the computation, with all weight points
-	double sum = 0.0;
-	double norm = 0.0;
-	double pi = 4.0*atan(1.0);
-	// No polydispersion supported in this model.
-	// Get the dispersion points for the radius
-	vector<WeightPoint> weights_radius;
-	radius.get_weights(weights_radius);
-	vector<WeightPoint> weights_edge_separation;
-	edge_separation.get_weights(weights_edge_separation);
-	// Perform the computation, with all weight points
-	double string_vol = 0.0;
-	double pearl_vol = 0.0;
-	double tot_vol = 0.0;
-	// Loop over core weight points
-	for(size_t i=0; i<weights_radius.size(); i++) {
-		dp.radius = weights_radius[i].value;
-		// Loop over thick_inter0 weight points
-		for(size_t j=0; j<weights_edge_separation.size(); j++) {
-			dp.edge_separation = weights_edge_separation[j].value;
-			pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
-			string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
-			tot_vol = (dp.num_pearls - 1.0) * string_vol;
-			tot_vol += dp.num_pearls * pearl_vol;
-			//Find  volume
-			// This may be a too much approximation
-			//Todo: decided whether or not we keep this calculation
-			sum += weights_radius[i].weight * weights_edge_separation[j].weight
-				* pow(3.0*tot_vol/4.0/pi,0.333333);
-			norm += weights_radius[i].weight * weights_edge_separation[j].weight;
-		}
-	}
-
-	if (norm != 0){
-		//return the averaged value
-		rad_out =  sum/norm;}
-	else{
-		//return normal value
-		pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
-		string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
-		tot_vol = (dp.num_pearls - 1.0) * string_vol;
-		tot_vol += dp.num_pearls * pearl_vol;
-
-		rad_out =  pow((3.0*tot_vol/4.0/pi), 0.33333);
-		}
-
-	return rad_out;
-
-}
+  PeralNecklaceParameters dp;
+
+  dp.scale = scale();
+  dp.radius = radius();
+  dp.edge_separation = edge_separation();
+  dp.thick_string = thick_string();
+  dp.num_pearls = num_pearls();
+
+  double rad_out = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double pi = 4.0*atan(1.0);
+  // No polydispersion supported in this model.
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+  vector<WeightPoint> weights_edge_separation;
+  edge_separation.get_weights(weights_edge_separation);
+  // Perform the computation, with all weight points
+  double string_vol = 0.0;
+  double pearl_vol = 0.0;
+  double tot_vol = 0.0;
+  // Loop over core weight points
+  for(size_t i=0; i<weights_radius.size(); i++) {
+    dp.radius = weights_radius[i].value;
+    // Loop over thick_inter0 weight points
+    for(size_t j=0; j<weights_edge_separation.size(); j++) {
+      dp.edge_separation = weights_edge_separation[j].value;
+      pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
+      string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
+      tot_vol = (dp.num_pearls - 1.0) * string_vol;
+      tot_vol += dp.num_pearls * pearl_vol;
+      //Find  volume
+      // This may be a too much approximation
+      //Todo: decided whether or not we keep this calculation
+      sum += weights_radius[i].weight * weights_edge_separation[j].weight
+          * pow(3.0*tot_vol/4.0/pi,0.333333);
+      norm += weights_radius[i].weight * weights_edge_separation[j].weight;
+    }
+  }
+
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
+    string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
+    tot_vol = (dp.num_pearls - 1.0) * string_vol;
+    tot_vol += dp.num_pearls * pearl_vol;
+
+    rad_out =  pow((3.0*tot_vol/4.0/pi), 0.33333);
+  }
+
+  return rad_out;
+
+}
Index: sansmodels/src/c_models/sld_cal.cpp
===================================================================
--- sansmodels/src/c_models/sld_cal.cpp	(revision 67424cdc019bde2f4cbc089828ff660e6f914c2c)
+++ sansmodels/src/c_models/sld_cal.cpp	(revision 0ba3b08c35684feff124c0657dad569d13717a18)
@@ -6,11 +6,45 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "sld_cal.h"
 
 extern "C" {
-	#include "sld_cal.h"
+#include "libmultifunc/librefl.h"
+}
+
+// Convenience structure
+typedef struct {
+    double fun_type;
+    double npts_inter;
+    double shell_num;
+    double nu_inter;
+    double sld_left;
+    double sld_right;
+} SLDCalParameters;
+
+/**
+ * Function to calculate sld
+ * @param pars: parameters
+ * @param x: independent param-value
+ * @return: sld value
+ */
+static double sld_cal_analytical_1D(SLDCalParameters *pars, double x) {
+  double fun, nsl, n_s, fun_coef, sld_l, sld_r, sld_out;
+  int fun_type;
+
+  fun = pars->fun_type;
+  nsl = pars->npts_inter;
+  n_s = pars->shell_num;
+  fun_coef = pars->nu_inter;
+  sld_l = pars-> sld_left;
+  sld_r = pars-> sld_right;
+
+  fun_type = floor(fun);
+
+  sld_out = intersldfunc(fun_type, nsl, n_s, fun_coef, sld_l, sld_r);
+
+  return sld_out;
 }
 
