ExhumeInterface/src/Integral.cpp

0001 //Simpson integrator, written to mirror CERNLIBS one
0002
0003 #include <cmath>
0004
0005 double dsimps_(double* F, double* A, double* B, int* N) {
0006   double a = *A;
0007   double b = *B;
0008   int n = *N;
0009   double deltaX = (b - a) / n;
0010   double I = F[0] + F[n];
0011
0012   for (int i = 1; i <= n / 2; i++) {
0013     int j = 2 * i;
0014     int k = 2 * i - 1;
0015     if (j < n)
0016       I += 2 * F[j];
0017     if (k < n)
0018       I += 4 * F[k];
0019   }
0020   return (I * deltaX / 3);
0021 }