/* With the functions in this file 1-dimensional definite integrals can be calculated using either - a Newton-Cotes algorithm : * trapezoidal rule * Simpson's rule - a Gauss-Legendre algorithm. */ #include #define MAXIT 20 /* this gives the max ( <= 32 !) number of iterations that is allowed */ #define FRACACC 1.0e-5 #define EPS 3.0e-11 /* relative precision */ #define FUNC(x) ((*func)(x)) /* func is a pointer to the function to be integrated */ #define FUNC2(x) ((*func2)(x)) /* func2 is a pointer to the function to be integrated */ #define FUNC3(x) ((*func3)(x)) /* func3 is a pointer to the function to be integrated */ #define FUNC4(x) ((*func4)(x)) /* func4 is a pointer to the function to be integrated */ #define FUNC5(x) ((*func5)(x)) /* func5 is a pointer to the function to be integrated */ #define FUNC6(x) ((*func6)(x)) /* func5 is a pointer to the function to be integrated */ #define FUNC7(x) ((*func7)(x)) /* func5 is a pointer to the function to be integrated */ /* The function int_trapzd() - which makes use of the function trapzd() - calculates the integral using the extended trapezoidal rule. See : Press et al., Numerical Recipes in C ,Cambridge 1992, 4.1 and 4.2 */ double int_trapzd(double (*func)(double x), double a, double b) { double previous = 1.0E-30; double integral = 1.0E-25; double fracacc; int j=1; fracacc = fabs( (previous - integral) / previous ); while( (fracacc > FRACACC) && (j <= MAXIT) ) { previous = integral; integral = trapzd(func, a, b, j); fracacc = fabs( (previous - integral) / previous ); j++; } return integral; } /* The function int_simp() - which makes use of the function trapzd() - calculates the integral using the extended simpson rule. See : Press et al., Numerical Recipes in C ,Cambridge 1992, 4.1 and 4.2 */ double int_simp(double (*func)(double x), double a, double b) { double previous = 1.0e-15; double integral = 1.0e-10; double integralterm1 , integralterm2 = 0.0; double fracacc; int j=1; fracacc = fabs( (previous - integral) / previous ); while( (fracacc > FRACACC) && (j <= MAXIT) ) { previous = integral; integralterm1 = trapzd(func, a, b, j++); integral = ((1.0 / 3.0) * (4.0 * integralterm1 - integralterm2)); integralterm2 = integralterm1; fracacc = fabs( (previous - integral) / previous ); } return integral; } double trapzd(double (*func)(double x), double a, double b, int n) { double coord, tnm, sum, del; static double s; long j,it; if (n == 1) { s = .5 * ( b - a) * (FUNC(a) + FUNC(b)); return s; } else { for (it = 1, j = 1; j < (n - 1); j++) it = it << 1; tnm = it; del = (b - a) / tnm; coord = a + .5 * del; for (sum=0.0,j=1; j<=it; j++,coord += del) sum += FUNC(coord); s = .5 * (s + (b - a) *sum / tnm); return s; } } /* The function int_gauss() calculates a one-dimensional definite integral using a gaussian quadrature formula (in the present case we use the GAUSS-LEGENDRE quadrature formula. - a and b are the lower and upper integration limits respectively, - n stands for the n-point quadrature formula (all polynomials up to degree (2 * n - 1) are integrated exactly), - x[1...n] and w[1...n] contain the abscissa and weights of the n-point Gauss-Legendre quadrature formula and are calculated by the function gauss_leg(). */ double int_gauss(double (*func)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x, *w; x = dvector(1, n); w = dvector(1, n); gauss_leg(a, b, x, w, n); for (j = 1; j <= n; j++) { sum += w[j] * FUNC(x[j]); } free_dvector(x, 1, n); free_dvector(w, 1, n); return sum; } double int_gauss2(double (*func2)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x2, *w2; x2 = dvector(1, n); w2 = dvector(1, n); gauss_leg(a, b, x2, w2, n); for (j = 1; j <= n; j++) { sum += w2[j] * FUNC2(x2[j]); } free_dvector(x2, 1, n); free_dvector(w2, 1, n); return sum; } double int_gauss3(double (*func3)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x3, *w3; x3 = dvector(1, n); w3 = dvector(1, n); gauss_leg(a, b, x3, w3, n); for (j = 1; j <= n; j++) { sum += w3[j] * FUNC3(x3[j]); } free_dvector(x3, 1, n); free_dvector(w3, 1, n); return sum; } double int_gauss4(double (*func4)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x4, *w4; x4 = dvector(1, n); w4 = dvector(1, n); gauss_leg(a, b, x4, w4, n); for (j = 1; j <= n; j++) { sum += w4[j] * FUNC4(x4[j]); } free_dvector(x4, 1, n); free_dvector(w4, 1, n); return sum; } double int_gauss5(double (*func5)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x5, *w5; x5 = dvector(1, n); w5 = dvector(1, n); gauss_leg(a, b, x5, w5, n); for (j = 1; j <= n; j++) { sum += w5[j] * FUNC5(x5[j]); } free_dvector(x5, 1, n); free_dvector(w5, 1, n); return sum; } double int_gauss6(double (*func6)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x6, *w6; x6 = dvector(1, n); w6 = dvector(1, n); gauss_leg(a, b, x6, w6, n); for (j = 1; j <= n; j++) { sum += w6[j] * FUNC6(x6[j]); } free_dvector(x6, 1, n); free_dvector(w6, 1, n); return sum; } double int_gauss7(double (*func7)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x7, *w7; x7 = dvector(1, n); w7 = dvector(1, n); gauss_leg(a, b, x7, w7, n); for (j = 1; j <= n; j++) { sum += w7[j] * FUNC7(x7[j]); } free_dvector(x7, 1, n); free_dvector(w7, 1, n); return sum; } /* The function gauss_leg() calculates the roots and weights for the GAUSS-LEGENDRE n-point quadrature formula. See : Press et al., Numerical Recipes in C, Cambridge UP 1992, section 4.5.*/ void gauss_leg(double x1, double x2, double *x, double *w, int n) { int m, j, i; double z1, z, xm, xl, pp, p3, p2, p1; m = (n + 1) / 2; xm = 0.5 * (x2 + x1); xl = 0.5 * (x2 - x1); for (i = 1; i <= m; i++) /* loop over desired roots */ { z = cos(3.141592654 * (i - 0.25)/(n + 0.5) ); /* Starting approximation to ith root */ do { p1 = 1.0; p2 = 0.0; for (j = 1; j <= n; j++) { p3 = p2; p2 = p1; p1 = ( (2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j ; } pp = n * (z * p1 - p2) / (z * z - 1.0); z1 = z; z = z1 - p1 / pp; /* Newton-Raphson interpolation formula */ } while (fabs(z - z1) > EPS); x[i] = xm - xl * z; /* Roots scaled to desired interval */ x[n + 1 - i] = xm + xl * z; /* symmetric counterparts of roots */ w[i] = 2.0 * xl / ( (1.0 - z * z) * pp * pp ); /* weights */ w[n + 1 - i] = w[i]; } } /* The function int_gauss_infinit() calculates a one-dimensional definite integral from zero to infinity using the GAUSS-LEGENDRE quadrature formula. - c is an arbitrary transformation factor that can be used to check the routine (the result should be independent of c), - n stands for the n-point quadrature formula (all polynomials up to degree (2 * n - 1) are integrated exactly), - x[1...n] and w[1...n] contain the abscissa and weights of the n-point Gauss-Legendre quadrature formula and are calculated by the function gauss_leg(). */ double int_gauss_infinit(double (*func)(double), double c, int n) { int j; double sum = 0.0, *q, *W; static double *x, *w; x = dvector(1, n); w = dvector(1, n); q = dvector(1, n); W = dvector(1, n); gauss_leg(- 0.99999, 0.99999, x, w, n); for (j = 1; j <= n; j++) { q[j] = c * tan(3.141592654 / 4.0 * (1.0 + x[j]) ); W[j] = c * 3.141592654 / 4.0 * w[j] * 1.0 / pow( cos(3.141592654 / 4.0 * (1.0 + x[j]) ), 2.0); sum += W[j] * FUNC(q[j]); } free_dvector(x, 1, n); free_dvector(w, 1, n); free_dvector(q, 1, n); free_dvector(W, 1, n); return sum; } /* The function gauss_leg_infinit() calculates the roots and weights for the GAUSS-LEGENDRE n-point quadrature formula for integration limits zero and infinity. */ void gauss_leg_infinit(double *q, double *W, double c, int n) { int j; static double *x, *w; x = dvector(1, n); w = dvector(1, n); gauss_leg(- 0.99999, 0.99999, x, w, n); for (j = 1; j <= n; j++) { q[j] = c * tan(3.141592654 / 4.0 * (1.0 + x[j]) ); W[j] = c * 3.141592654 / 4.0 * w[j] * 1.0 / pow( cos(3.141592654 / 4.0 * (1.0 + x[j]) ), 2.0); } free_dvector(x, 1, n); free_dvector(w, 1, n); }