Index: sasmodels/models/pringle.c
===================================================================
--- sasmodels/models/pringle.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/pringle.c	(revision c047acf8fad567cad9bc9c45d3e56d4805a8733f)
@@ -1,4 +1,3 @@
-double form_volume(double radius,
-          double thickness);
+double form_volume(double radius, double thickness, double alpha, double beta);
 
 double Iq(double q,
@@ -10,171 +9,123 @@
           double sld_solvent);
 
+
 static
-double pringleC(double radius,
-                double alpha,
-                double beta,
-                double q,
-                double phi,
-                double n) {
-
-    double va, vb;
-    double bessargs, cosarg, bessargcb;
-    double r, retval, yyy;
-
-
-    va = 0;
-    vb = radius;
+void _integrate_bessel(
+    double radius,
+    double alpha,
+    double beta,
+    double q_sin_psi,
+    double q_cos_psi,
+    double n,
+    double *Sn,
+    double *Cn)
+{
+    // translate gauss point z in [-1,1] to a point in [0, radius]
+    const double zm = 0.5*radius;
+    const double zb = 0.5*radius;
 
     // evaluate at Gauss points
-    // remember to index from 0,size-1
+    double sumS = 0.0;		// initialize integral
+    double sumC = 0.0;		// initialize integral
+    double r;
+    for (int i=0; i < 76; i++) {
+        r = Gauss76Z[i]*zm + zb;
 
-    double summ = 0.0;		// initialize integral
-    int ii = 0;
-    do {
-        // Using 76 Gauss points
-        r = (Gauss76Z[ii] * (vb - va) + vb + va) / 2.0;
+        const double qrs = r*q_sin_psi;
+        const double qrrc = r*r*q_cos_psi;
 
-        bessargs = q*r*sin(phi);
-        cosarg = q*r*r*alpha*cos(phi);
-        bessargcb = q*r*r*beta*cos(phi);
+        double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);
+        double S, C;
+        SINCOS(alpha*qrrc, S, C);
+        sumS += y*S;
+        sumC += y*C;
+    }
 
-        yyy = Gauss76Wt[ii]*r*cos(cosarg)
-                *sas_JN(n, bessargcb)
-                *sas_JN(2*n, bessargs);
-        summ += yyy;
-
-        ii += 1;
-    } while (ii < N_POINTS_76);			// 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;
+    #if 1
+    // TODO: should the normalization be to R^2 or the last value, r^2
+    // the gaussian window does not go all the way from 0 to 1.
+    //radius = Gauss76Z[75] * zm + zb;
+    *Sn = zm*sumS / (r*r);
+    *Cn = zm*sumC / (r*r);
+    #else
+    *Sn = zm*sumS / (radius*radius);
+    *Cn = zm*sumC / (radius*radius);
+    #endif
 }
 
 static
-double pringleS(double radius,
-                double alpha,
-                double beta,
-                double q,
-                double phi,
-                double n) {
-
-    double va, vb, summ;
-    double bessargs, sinarg, bessargcb;
-    double r, retval, yyy;
-    // set up the integration
-    // end points and weights
-
-    va = 0;
-    vb = radius;
-
-    // evaluate at Gauss points
-    // remember to index from 0,size-1
-
-    summ = 0.0;		// initialize integral
-    int 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*alpha*cos(phi);
-        bessargcb = q*r*r*beta*cos(phi);
-
-        yyy = Gauss76Wt[ii]*r*sin(sinarg)
-                    *sas_JN(n, bessargcb)
-                    *sas_JN(2*n, bessargs);
-        summ += yyy;
-
-        ii += 1;
-    } while (ii < N_POINTS_76);
-
-    // 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;
+double _sum_bessel_orders(
+    double radius,
+    double alpha,
+    double beta,
+    double q_sin_psi,
+    double q_cos_psi)
+{
+    //calculate sum term from n = -3 to 3
+    //Note 1:
+    //    S_n(-x) = (-1)^S_n(x)
+    //    => S_n^2(-x) = S_n^2(x),
+    //    => sum_-k^k Sk = S_0^2 + 2*sum_1^kSk^2
+    //Note 2:
+    //    better precision to sum terms from smaller to larger
+    //    though it doesn't seem to make a difference in this case.
+    double Sn, Cn, sum;
+    sum = 0.0;
+    for (int n=3; n>0; n--) {
+      _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, n, &Sn, &Cn);
+      sum += 2.0*(Sn*Sn + Cn*Cn);
+    }
+    _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, 0, &Sn, &Cn);
+    sum += Sn*Sn+ Cn*Cn;
+    return sum;
 }
 
 static
-double _kernel(double thickness,
-               double radius,
-               double alpha,
-               double beta,
-               double q,
-               double phi) {
+double _integrate_psi(
+    double q,
+    double radius,
+    double thickness,
+    double alpha,
+    double beta)
+{
+    // translate gauss point z in [-1,1] to a point in [0, pi/2]
+    const double zm = M_PI_4;
+    const double zb = M_PI_4;
 
-    const double sincarg = q * thickness * cos(phi) / 2.0;
-    const double sincterm = pow(sin(sincarg) / sincarg, 2.0);
+    double sum = 0.0;
+    for (int i = 0; i < 76; i++) {
+        double psi = Gauss76Z[i]*zm + zb;
+        double sin_psi, cos_psi;
+        SINCOS(psi, sin_psi, cos_psi);
+        double bessel_term = _sum_bessel_orders(radius, alpha, beta, q*sin_psi, q*cos_psi);
+        double sinc_term = square(sinc(q * thickness * cos_psi / 2.0));
+        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term;
+        sum += Gauss76Wt[i] * pringle_kernel;
+    }
 
-    //calculate sum term from n = -3 to 3
-    double sumterm = 0.0;
-    for (int nn = -3; nn <= 3; nn++) {
-        double powc = pringleC(radius, alpha, beta, q, phi, nn);
-        double pows = pringleS(radius, alpha, beta, q, phi, nn);
-        sumterm += pow(powc, 2.0) + pow(pows, 2.0);
-    }
-    double retval = 4.0 * sin(phi) * sumterm * sincterm;
-
-    return retval;
-
+    return zm * sum;
 }
 
-static double pringles_kernel(double q,
-          double radius,
-          double thickness,
-          double alpha,
-          double beta,
-          double sld_pringle,
-          double sld_solvent)
+double form_volume(double radius, double thickness, double alpha, double beta)
 {
-
-    //upper and lower integration limits
-    const double lolim = 0.0;
-    const double uplim = M_PI / 2.0;
-
-    double summ = 0.0;			//initialize integral
-
-    double delrho = sld_pringle - sld_solvent; //make contrast term
-
-    for (int i = 0; i < N_POINTS_76; i++) {
-        double phi = (Gauss76Z[i] * (uplim - lolim) + uplim + lolim) / 2.0;
-        summ += Gauss76Wt[i] * _kernel(thickness, radius, alpha, beta, q, phi);
-    }
-
-    double answer = (uplim - lolim) / 2.0 * summ;
-    answer *= delrho*delrho;
-
-    return answer;
+    // TODO: Normalize by form volume
+    //return M_PI*radius*radius*thickness;
+    return 1.0;
 }
 
-double form_volume(double radius,
-        double thickness){
-
-        return 1.0;
+double Iq(
+    double q,
+    double radius,
+    double thickness,
+    double alpha,
+    double beta,
+    double sld_pringle,
+    double sld_solvent)
+{
+    double form = _integrate_psi(q, radius, thickness, alpha, beta);
+    double contrast = sld_pringle - sld_solvent;
+    double volume = M_PI*radius*radius*thickness;
+    // TODO: If normalize by form volume, need an extra volume here
+    //return 1.0e-4*form * square(contrast * volume);
+    return 1.0e-4*form * square(contrast) * volume;
 }
-
-double Iq(double q,
-          double radius,
-          double thickness,
-          double alpha,
-          double beta,
-          double sld_pringle,
-          double sld_solvent)
-{
-    const double form = pringles_kernel(q,
-                  radius,
-                  thickness,
-                  alpha,
-                  beta,
-                  sld_pringle,
-                  sld_solvent);
-
-    return 1.0e-4*form*M_PI*radius*radius*thickness;
-}
Index: sasmodels/models/pringle.py
===================================================================
--- sasmodels/models/pringle.py	(revision 42356c8890c04b71169941b12a36fc0e07e2034b)
+++ sasmodels/models/pringle.py	(revision c047acf8fad567cad9bc9c45d3e56d4805a8733f)
@@ -66,6 +66,6 @@
     ["radius",      "Ang",         60.0,   [0, inf],    "volume", "Pringle radius"],
     ["thickness",   "Ang",         10.0,   [0, inf],    "volume", "Thickness of pringle"],
-    ["alpha",       "",            0.001,  [-inf, inf], "", "Curvature parameter alpha"],
-    ["beta",        "",            0.02,   [-inf, inf], "", "Curvature paramter beta"],
+    ["alpha",       "",            0.001,  [-inf, inf], "volume", "Curvature parameter alpha"],
+    ["beta",        "",            0.02,   [-inf, inf], "volume", "Curvature paramter beta"],
     ["sld_pringle", "1e-6/Ang^2",  1.0,    [-inf, inf], "sld", "Pringle sld"],
     ["sld_solvent", "1e-6/Ang^2",  6.3,    [-inf, inf], "sld", "Solvent sld"]
@@ -73,8 +73,9 @@
 # pylint: enable=bad-whitespace, line-too-long
 
+
 source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", \
           "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"]
 
-def ER(radius, thickness):
+def ER(radius, thickness, alpha, beta):
     """
         Return equivalent radius (ER)
