Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision 6831fa0e3e5df0adb39649e7ec1ac9ce83f8c975)
+++ sasmodels/models/stacked_disks.c	(revision 3ac4e1bbf80dc45fabe12c92c0d88f895d96152e)
@@ -27,5 +27,5 @@
 
 static
-double _kernel(double qq,
+double _kernel(double q,
                double radius,
                double core_sld,
@@ -41,44 +41,43 @@
 
 {
-	// qq is the q-value for the calculation (1/A)
-	// radius is the core radius of the cylinder (A)
-	// *_sld are the respective SLD's
-	// halfheight is the *Half* CORE-LENGTH of the cylinder = L (A)
-	// zi is the dummy variable for the integration (x in Feigin's notation)
+    // q is the q-value for the calculation (1/A)
+    // radius is the core radius of the cylinder (A)
+    // *_sld are the respective SLD's
+    // halfheight is the *Half* CORE-LENGTH of the cylinder = L (A)
+    // zi is the dummy variable for the integration (x in Feigin's notation)
 
-	const double besarg1 = qq*radius*sin_alpha;
-	//const double besarg2 = qq*radius*sin_alpha;
+    const double besarg1 = q*radius*sin_alpha;
+    //const double besarg2 = q*radius*sin_alpha;
 
-	const double sinarg1 = qq*halfheight*cos_alpha;
-	const double sinarg2 = qq*(halfheight+thick_layer)*cos_alpha;
+    const double sinarg1 = q*halfheight*cos_alpha;
+    const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha;
 
-	const double be1 = sas_J1c(besarg1);
-	//const double be2 = sas_J1c(besarg2);
-	const double be2 = be1;
-	const double si1 = sinc(sinarg1);
-	const double si2 = sinc(sinarg2);
+    const double be1 = sas_J1c(besarg1);
+    //const double be2 = sas_J1c(besarg2);
+    const double be2 = be1;
+    const double si1 = sinc(sinarg1);
+    const double si2 = sinc(sinarg2);
 
-	const double dr1 = (core_sld-solvent_sld);
-	const double dr2 = (layer_sld-solvent_sld);
-	const double area = M_PI*radius*radius;
-	const double totald = 2.0*(thick_layer+halfheight);
+    const double dr1 = core_sld - solvent_sld;
+    const double dr2 = layer_sld - solvent_sld;
+    const double area = M_PI*radius*radius;
+    const double totald = 2.0*(thick_layer + halfheight);
 
-	const double t1 = area*(2.0*halfheight)*dr1*(si1)*(be1);
-	const double t2 = area*dr2*(totald*si2-2.0*halfheight*si1)*(be2);
+    const double t1 = area * (2.0*halfheight) * dr1 * si1 * be1;
+    const double t2 = area * dr2 * (totald*si2 - 2.0*halfheight*si1) * be2;
 
+    double pq = square(t1 + t2);
 
-	double retval =((t1+t2)*(t1+t2));
+    // loop for the structure factor S(q)
+    double qd_cos_alpha = q*d*cos_alpha;
+    double debye_arg = -0.5*square(qd_cos_alpha*sigma_dnn);
+    double sq=0.0;
+    for (int kk=1; kk<n_stacking; kk++) {
+        sq += (n_stacking-kk) * cos(qd_cos_alpha*kk) * exp(debye_arg*kk);
+    }
+    // end of loop for S(q)
+    sq = 1.0 + 2.0*sq/n_stacking;
 
-	// loop for the structure facture S(q)
-	double sqq=0.0;
-	for(int kk=1;kk<n_stacking;kk+=1) {
-		double qd_cos_alpha = qq*d*cos_alpha;
-		double dexpt=square(qd_cos_alpha*sigma_dnn)*kk/2.0;
-		sqq += (n_stacking-kk)*cos(qd_cos_alpha*kk)*exp(-1.*dexpt);
-	}
-	// end of loop for S(q)
-	sqq=1.0+2.0*sqq/n_stacking;
-
-	return retval * sqq;
+    return pq * sq;
 }
 
@@ -95,37 +94,35 @@
                             double solvent_sld)
 {
-/*	StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
+/*    StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
 like clay platelets that are not exfoliated
 */
-	double summ = 0.0;	//initialize integral
+    double summ = 0.0;    //initialize integral
 
-	double d = 2.0*thick_layer+thick_core;
-	double halfheight = 0.5*thick_core;
+    double d = 2.0*thick_layer+thick_core;
+    double halfheight = 0.5*thick_core;
 
-	for(int i=0; i<N_POINTS_76; i++) {
-		double zi = (Gauss76Z[i] + 1.0)*M_PI_4;
-		double sin_alpha, cos_alpha; // slots to hold sincos function output
-		SINCOS(zi, sin_alpha, cos_alpha);
-		double yyy = _kernel(q,
-		                   radius,
-		                   core_sld,
-		                   layer_sld,
-		                   solvent_sld,
-		                   halfheight,
-		                   thick_layer,
-		                   sin_alpha,
-		                   cos_alpha,
-		                   sigma_dnn,
-		                   d,
-		                   n_stacking);
-		summ += Gauss76Wt[i] * yyy * sin_alpha;
-	}
+    for(int i=0; i<N_POINTS_76; i++) {
+        double zi = (Gauss76Z[i] + 1.0)*M_PI_4;
+        double sin_alpha, cos_alpha; // slots to hold sincos function output
+        SINCOS(zi, sin_alpha, cos_alpha);
+        double yyy = _kernel(q,
+                           radius,
+                           core_sld,
+                           layer_sld,
+                           solvent_sld,
+                           halfheight,
+                           thick_layer,
+                           sin_alpha,
+                           cos_alpha,
+                           sigma_dnn,
+                           d,
+                           n_stacking);
+        summ += Gauss76Wt[i] * yyy * sin_alpha;
+    }
 
-	double answer = M_PI_4*summ;
+    double answer = M_PI_4*summ;
 
-	//Convert to [cm-1]
-	answer *= 1.0e-4;
-
-	return answer;
+    //Convert to [cm-1]
+    return 1.0e-4*answer;
 }
 
