Index: sansmodels/src/c_models/pringles.cpp
===================================================================
--- sansmodels/src/c_models/pringles.cpp	(revision a9fec15f74f829e85ab0ebd3b7f671ad70e3c00a)
+++ sansmodels/src/c_models/pringles.cpp	(revision a9fec15f74f829e85ab0ebd3b7f671ad70e3c00a)
@@ -0,0 +1,350 @@
+/**
+ * PringlesModel
+ *
+ * (c) 2013 / Andrew J Jackson / andrew.jackson@esss.se
+ *
+ * Model for "pringles" particle from K. Edler @ Bath University
+ *
+ */
+#include <math.h>
+#include "parameters.hh"
+#include "pringles.h"
+using namespace std;
+
+extern "C"
+{
+#include "GaussWeights.h"
+#include "libStructureFactor.h"
+}
+
+// Convenience parameter structure
+typedef struct {
+    double scale;
+    double radius;
+    double thickness;
+    double alpha;
+    double beta;
+    double sldCyl;
+    double sldSolv;
+    double background;
+    double cyl_theta;
+    double cyl_phi;
+} PringleParameters;
+
+
+PringlesModel::PringlesModel() {
+	scale = Parameter(1.0);
+	radius = Parameter(60.0,true);
+	radius.set_min(0.0);
+	thickness = Parameter(10.0,true);
+	thickness.set_min(0.0);
+	alpha = Parameter(0.001,true);
+	beta = Parameter(0.02,true);
+	sld_pringle = Parameter(1.0e-6);
+	sld_solvent = Parameter(6.35e-6);
+	background = Parameter(0.0);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * @param q: q-value
+ * @return: function value
+ */
+double PringlesModel::operator()(double q) {
+	double dp[8];
+	// Fill parameter array
+	// Add the background after averaging
+	dp[0] = scale();
+	dp[1] = radius();
+	dp[2] = thickness();
+	dp[3] = alpha();
+	dp[4] = beta();
+	dp[5] = sld_pringle();
+	dp[6] = sld_solvent();
+	dp[7] = 0.0;
+
+	// Get the dispersion points for the radius
+	vector<WeightPoint> weights_rad;
+	radius.get_weights(weights_rad);
+
+	// Get the dispersion points for the thickness
+	vector<WeightPoint> weights_thick;
+	thickness.get_weights(weights_thick);
+
+	// Get the dispersion points for alpha
+	vector<WeightPoint> weights_alpha;
+	alpha.get_weights(weights_alpha);
+
+	// Get the dispersion points for beta
+	vector<WeightPoint> weights_beta;
+	beta.get_weights(weights_beta);
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+	double volnorm = 0.0;
+	double vol = 0.0;
+	double Pi;
+
+	Pi = 4.0 * atan(1.0);
+
+	// Loop over alpha weight points
+	for (size_t i = 0; i < weights_alpha.size(); i++) {
+		dp[3] = weights_alpha[i].value;
+
+		//Loop over beta weight points
+		for (size_t j = 0; j < weights_beta.size(); j++) {
+			dp[4] = weights_beta[j].value;
+
+			// Loop over thickness weight points
+			for (size_t k = 0; k < weights_thick.size(); k++) {
+				dp[2] = weights_thick[k].value;
+
+				// Loop over radius weight points
+				for (size_t l = 0; l < weights_rad.size(); l++) {
+					dp[1] = weights_rad[l].value;
+					sum += weights_rad[l].weight * weights_thick[k].weight
+							* weights_alpha[i].weight * weights_beta[j].weight
+							* pringle_form(dp, q);
+					//Find average volume
+					vol += weights_rad[l].weight * weights_thick[k].weight * Pi * pow(weights_rad[l].value, 2) * weights_thick[k].value;
+					volnorm += weights_rad[l].weight * weights_thick[k].weight;
+					norm +=  weights_rad[l].weight * weights_thick[k].weight * weights_alpha[i].weight * weights_beta[j].weight;
+
+				}
+			}
+		}
+	}
+
+	if (vol > 0.0 && norm > 0.0) {
+		//normalize by avg volume
+		sum = sum * (vol/volnorm);
+		return sum/norm + background();
+	} else {
+		return 0.0;
+	}
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double PringlesModel::operator()(double qx, double qy) {
+	double q = sqrt(qx * qx + qy * qy);
+	return (*this).operator()(q);
+}
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the model
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double PringlesModel::evaluate_rphi(double q, double phi) {
+	return (*this).operator()(q);
+}
+
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double PringlesModel::calculate_ER() {
+	PringleParameters dp;
+
+	dp.radius     = radius();
+	dp.thickness  = thickness();
+	double rad_out = 0.0;
+
+	// Perform the computation, with all weight points
+	double sum = 0.0;
+	double norm = 0.0;
+
+	// Get the dispersion points for the major shell
+	vector<WeightPoint> weights_thick;
+	thickness.get_weights(weights_thick);
+
+	// Get the dispersion points for the minor shell
+	vector<WeightPoint> weights_radius ;
+	radius.get_weights(weights_radius);
+
+	// Loop over major shell weight points
+	for(int i=0; i< (int)weights_thick.size(); i++) {
+		dp.thickness = weights_thick[i].value;
+		for(int k=0; k< (int)weights_radius.size(); k++) {
+			dp.radius = weights_radius[k].value;
+			//Note: output of "DiamCyl(dp.thick,dp.radius)" is DIAMETER.
+			sum +=weights_thick[i].weight
+				* weights_radius[k].weight*DiamCyl(dp.thickness,dp.radius)/2.0;
+			norm += weights_thick[i].weight* weights_radius[k].weight;
+		}
+	}
+	if (norm != 0){
+		//return the averaged value
+		rad_out =  sum/norm;}
+	else{
+		//return normal value
+		//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+		rad_out = DiamCyl(dp.thickness,dp.radius)/2.0;}
+
+	return rad_out;
+}
+/**
+ * Function to calculate particle volume/total volume for shape models:
+ *    Most case returns 1 but for example for the vesicle model it is
+ *    (total volume - core volume)/total volume
+ *    (< 1 depending on the thickness).
+ * @return: effective radius value
+ */
+double PringlesModel::calculate_VR() {
+	return 1.0;
+}
+
+/*
+ * Useful work functions start here!
+ *
+ */
+
+static double pringle_form(double dp[], double q) {
+
+	double Pi;
+	int nord = 76, i=0;			//order of integration
+	double uplim, lolim;		//upper and lower integration limits
+	double summ, phi, yyy, answer, vcyl;			//running tally of integration
+	double delrho;
+
+	Pi = 4.0 * atan(1.0);
+	lolim = 0.0;
+	uplim = Pi / 2.0;
+
+	summ = 0.0;			//initialize integral
+
+	delrho = dp[5] - dp[6] ; //make contrast term
+
+	for (i = 0; i < nord; i++) {
+		phi = (Gauss76Z[i] * (uplim - lolim) + uplim + lolim) / 2.0;
+		yyy = Gauss76Wt[i] * pringle_kernel(dp, q, phi);
+		summ += yyy;
+	}
+
+	answer = (uplim - lolim) / 2.0 * summ;
+
+	answer *= delrho*delrho;
+
+    //normalize by cylinder volume
+	//vcyl=Pi*dp[1]*dp[1]*dp[2];
+	//answer *= vcyl;
+
+    //convert to [cm-1]
+	answer *= 1.0e8;
+
+    //Scale by volume fraction
+	answer *= dp[0];
+
+	return answer;
+}
+
+static double pringle_kernel(double dp[], double q, double phi) {
+
+	double sumterm, sincarg, sincterm, nn, retval;
+
+	sincarg = q * dp[2] * cos(phi) / 2.0; //dp[2] = thickness
+	sincterm = pow(sin(sincarg) / sincarg, 2.0);
+
+	//calculate sum term from n = -3 to 3
+	sumterm = 0;
+	for (nn = -3; nn <= 3; nn = nn + 1) {
+		sumterm = sumterm
+				+ (pow(pringleC(dp, q, phi, nn), 2.0)
+						+ pow(pringleS(dp, q, phi, nn), 2.0));
+	}
+
+	retval = 4.0 * sin(phi) * sumterm * sincterm;
+
+	return retval;
+
+}
+
+static double pringleC(double dp[], double q, double phi, double n) {
+
+	double nord, va, vb, summ;
+	double bessargs, cosarg, bessargcb;
+	double r, retval, yyy;
+	int ii;
+	// set up the integration
+	// end points and weights
+	nord = 76;
+	va = 0;
+	vb = dp[1]; //radius
+
+	// evaluate at Gauss points
+	// remember to index from 0,size-1
+
+	summ = 0.0;		// initialize integral
+	ii = 0;
+	do {
+		// Using 76 Gauss points
+		r = (Gauss76Z[ii] * (vb - va) + vb + va) / 2.0;
+
+		bessargs = q * r * sin(phi);
+		cosarg = q * r * r * dp[3] * cos(phi);
+		bessargcb = q * r * r * dp[4] * cos(phi);
+
+		yyy = Gauss76Wt[ii] * r * cos(cosarg) * jn(n, bessargcb)
+				* jn(2 * n, bessargs);
+		summ += yyy;
+
+		ii += 1;
+	} while (ii < nord);			// end of loop over quadrature points
+	//
+	// calculate value of integral to return
+
+	retval = (vb - va) / 2.0 * summ;
+
+	retval = retval / pow(r, 2.0);
+
+	return retval;
+}
+
+static double pringleS(double dp[], double q, double phi, double n) {
+
+	double nord, va, vb, summ;
+	double bessargs, sinarg, bessargcb;
+	double r, retval, yyy;
+	int ii;
+	// set up the integration
+	// end points and weights
+	nord = 76;
+	va = 0;
+	vb = dp[1]; //radius
+
+	// evaluate at Gauss points
+	// remember to index from 0,size-1
+
+	summ = 0.0;		// initialize integral
+	ii = 0;
+	do {
+		// Using 76 Gauss points
+		r = (Gauss76Z[ii] * (vb - va) + vb + va) / 2.0;
+
+		bessargs = q * r * sin(phi);
+		sinarg = q * r * r * dp[3] * cos(phi);
+		bessargcb = q * r * r * dp[4] * cos(phi);
+
+		yyy = Gauss76Wt[ii] * r * sin(sinarg) * jn(n, bessargcb)
+				* jn(2 * n, bessargs);
+
+		summ += yyy;
+
+		ii += 1;
+	} while (ii < nord);			// end of loop over quadrature points
+	//
+	// calculate value of integral to return
+
+	retval = (vb - va) / 2.0 * summ;
+
+	retval = retval / pow(r, 2.0);
+
+	return retval;
+}
