Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
+++ sasmodels/models/stacked_disks.c	(revision 6831fa0e3e5df0adb39649e7ec1ac9ce83f8c975)
@@ -14,4 +14,16 @@
           double solvent_sld);
 
+double Iqxy(double qx, double qy,
+          double thick_core,
+          double thick_layer,
+          double radius,
+          double n_stacking,
+          double sigma_dnn,
+          double core_sld,
+          double layer_sld,
+          double solvent_sld,
+          double theta,
+          double phi);
+
 static
 double _kernel(double qq,
@@ -22,5 +34,6 @@
                double halfheight,
                double thick_layer,
-               double zi,
+               double sin_alpha,
+               double cos_alpha,
                double sigma_dnn,
                double d,
@@ -34,19 +47,20 @@
 	// zi is the dummy variable for the integration (x in Feigin's notation)
 
-	const double besarg1 = qq*radius*sin(zi);
-	const double besarg2 = qq*radius*sin(zi);
+	const double besarg1 = qq*radius*sin_alpha;
+	//const double besarg2 = qq*radius*sin_alpha;
 
-	const double sinarg1 = qq*halfheight*cos(zi);
-	const double sinarg2 = qq*(halfheight+thick_layer)*cos(zi);
+	const double sinarg1 = qq*halfheight*cos_alpha;
+	const double sinarg2 = qq*(halfheight+thick_layer)*cos_alpha;
 
 	const double be1 = sas_J1c(besarg1);
-	const double be2 = sas_J1c(besarg2);
-	const double si1 = sin(sinarg1)/sinarg1;
-	const double si2 = sin(sinarg2)/sinarg2;
+	//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 totald = 2.0*(thick_layer+halfheight);
 
 	const double t1 = area*(2.0*halfheight)*dr1*(si1)*(be1);
@@ -54,19 +68,17 @@
 
 
-	double retval =((t1+t2)*(t1+t2))*sin(zi);
+	double retval =((t1+t2)*(t1+t2));
 
 	// loop for the structure facture S(q)
 	double sqq=0.0;
 	for(int kk=1;kk<n_stacking;kk+=1) {
-		double dexpt=qq*cos(zi)*qq*cos(zi)*d*d*sigma_dnn*sigma_dnn*kk/2.0;
-		sqq=sqq+(n_stacking-kk)*cos(qq*cos(zi)*d*kk)*exp(-1.*dexpt);
+		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;
 
-	retval *= sqq;
-
-	return(retval);
+	return retval * sqq;
 }
 
@@ -88,11 +100,12 @@
 	double summ = 0.0;	//initialize integral
 
-	double d=2.0*thick_layer+thick_core;
-	double halfheight = thick_core/2.0;
+	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.0;
-		double yyy = Gauss76Wt[i] *
-                    _kernel(q,
+	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,
@@ -101,12 +114,13 @@
 		                   halfheight,
 		                   thick_layer,
-		                   zi,
+		                   sin_alpha,
+		                   cos_alpha,
 		                   sigma_dnn,
 		                   d,
 		                   n_stacking);
-		summ += yyy;
+		summ += Gauss76Wt[i] * yyy * sin_alpha;
 	}
 
-	double answer = M_PI/4.0*summ;
+	double answer = M_PI_4*summ;
 
 	//Convert to [cm-1]
@@ -116,69 +130,10 @@
 }
 
-static double stacked_disks_kernel_2d(double q, double q_x, double q_y,
-                            double thick_core,
-                            double thick_layer,
-                            double radius,
-                            double n_stacking,
-                            double sigma_dnn,
-                            double core_sld,
-                            double layer_sld,
-                            double solvent_sld,
-                            double theta,
-                            double phi)
-{
-
-    double ct, st, cp, sp;
-
-    //convert angle degree to radian
-    theta = theta * M_PI/180.0;
-    phi = phi * M_PI/180.0;
-
-    SINCOS(theta, st, ct);
-    SINCOS(phi, sp, cp);
-
-    // silence compiler warnings about unused variable
-    (void) sp;
-
-    // parallelepiped orientation
-    const double cyl_x = st * cp;
-    const double cyl_y = st * sp;
-
-    // Compute the angle btw vector q and the
-    // axis of the parallelepiped
-    const double cos_val = cyl_x*q_x + cyl_y*q_y;
-
-    // Note: cos(alpha) = 0 and 1 will get an
-    // undefined value from Stackdisc_kern
-    double alpha = acos( cos_val );
-
-    // Call the IGOR library function to get the kernel
-    double d = 2 * thick_layer + thick_core;
-    double halfheight = thick_core/2.0;
-    double answer = _kernel(q,
-                     radius,
-                     core_sld,
-                     layer_sld,
-                     solvent_sld,
-                     halfheight,
-                     thick_layer,
-                     alpha,
-                     sigma_dnn,
-                     d,
-                     n_stacking);
-
-    answer /= sin(alpha);
-    //convert to [cm-1]
-    answer *= 1.0e-4;
-
-    return answer;
-}
-
 double form_volume(double thick_core,
                    double thick_layer,
                    double radius,
                    double n_stacking){
-    double d = 2 * thick_layer + thick_core;
-    return acos(-1.0) * radius * radius * d * n_stacking;
+    double d = 2.0 * thick_layer + thick_core;
+    return M_PI * radius * radius * d * n_stacking;
 }
 
@@ -203,2 +158,41 @@
                     solvent_sld);
 }
+
+
+double
+Iqxy(double qx, double qy,
+     double thick_core,
+     double thick_layer,
+     double radius,
+     double n_stacking,
+     double sigma_dnn,
+     double core_sld,
+     double layer_sld,
+     double solvent_sld,
+     double theta,
+     double phi)
+{
+    double q, sin_alpha, cos_alpha;
+    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+
+    double d = 2.0 * thick_layer + thick_core;
+    double halfheight = 0.5*thick_core;
+    double answer = _kernel(q,
+                     radius,
+                     core_sld,
+                     layer_sld,
+                     solvent_sld,
+                     halfheight,
+                     thick_layer,
+                     sin_alpha,
+                     cos_alpha,
+                     sigma_dnn,
+                     d,
+                     n_stacking);
+
+    //convert to [cm-1]
+    answer *= 1.0e-4;
+
+    return answer;
+}
+
