diff --git a/Qk.h b/Qk.h
index 7f392f1..326a342 100644
--- a/Qk.h
+++ b/Qk.h
@@ -40,61 +40,43 @@ double * QK(double Energy_keV, double radius_cm, double distance_cm, double thic
   double delx2 = (gamma-alpha)/div;
   
   for(int i = 0; i<=div; i++){
-
-    int J = i % 2;
-    if( J == 0 ){
-      A = 2.0;
-    }else{
-      A = 4.0;
-      beta1 = i * delx1;
-      beta2 = alpha + i * delx2;
-    }
     
     if( i == 0 || i == div ) {
       A = 1.0;
       beta1 = i * delx1;
       beta2 = alpha + i * delx2;
+    }else{
+      A = ( i % 2 == 0 ? 2.0 : 4.0 );
+      if( i % 2 == 1 ){
+        beta1 = i * delx1;
+        beta2 = alpha + i * delx2;
+      }
     }
 
-    double cosb = cos(beta1);
-    double sinb = sin(beta1);
-    double ex1 = exp( - 1.0 * Tau * thickness_cm / cosb );
-    double term3 = (1-ex1) * sinb * A * delx1;
-    double term1 = LegendreP(2, beta1) * term3;
-    double term2 = LegendreP(4, beta1) * term3;
+    double ex1 = exp( - 1.0 * Tau * thickness_cm / cos(beta1) );
+    double term1 = A * (1-ex1) * sin(beta1) * delx1;
 
-    sum1 += term1;
-    sum2 += term2;
-    sum3 += term3;
+    sum1 += LegendreP(2, beta1) * term1;
+    sum2 += LegendreP(4, beta1) * term1;
+    sum3 += term1;
 
-    cosb = cos(beta2);
-    sinb = sin(beta2);
-    double ex2 = exp( -1 * Tau * (radius_cm / sinb - distance_cm /cosb) );
-    double term6 = A * (1-ex2) * sinb * delx2;
-    double term4 = LegendreP(2, beta2) * term6;
-    double term5 = LegendreP(4, beta2) * term6;
+    double ex2 = exp( -1 * Tau * (radius_cm / sin(beta2) - distance_cm /cos(beta2)) );
+    double term2 = A * (1-ex2) * sin(beta2) * delx2;
 
-    sum4 = sum4 +term4;
-    sum5 = sum5 +term5;
-    sum6 = sum6 +term6;
+    sum4 += LegendreP(2, beta2) * term2;
+    sum5 += LegendreP(4, beta2) * term2;
+    sum6 += term2;
 
-    //if( i % 75 == 0) printf("%4d | %10.6f, %10.6f %d | %10.6f, %10.6f %10.6f \n", i, beta, A, J, sum1, sum2, sum3);
+    //if( i % 75 == 0) printf("%4d | %10.6f, %10.6f, %10.6f %d | %10.6f, %10.6f %10.6f %10.6f, %10.6f %10.6f \n", i, beta1, beta2, A, J, sum1, sum2, sum3, sum4, sum5, sum6);
   }
-    
-  double ans1 = sum1/3;
-  double ans2 = sum2/3;
-  double ans3 = sum3/3;  
-  double ans4 = sum4/3;
-  double ans5 = sum5/3;
-  double ans6 = sum6/3;
-
+  
   double * Qk = new double[2];
-  Qk[0]  = (ans1+ans4)/(ans3+ans6);
-  Qk[1]  = (ans2+ans5)/(ans3+ans6);
+  Qk[0]  = (sum1+sum4)/(sum3+sum6);
+  Qk[1]  = (sum2+sum5)/(sum3+sum6);
   
   printf("--------------\n");
-  printf("  QD2 = %lf\n",Qk[0]);
-  printf("  QD4 = %lf\n",Qk[1]);
+  printf("  QD2 = %.20f\n",Qk[0]);
+  printf("  QD4 = %.20f\n",Qk[1]);
   printf("--------------\n");
 
   return Qk;