9 * Compute the Nth stage of refinement of an extended trapezoidal
10 * rule. FUNC is input as a pointer to a function to be integrated
11 * between limits A and B. When called with N = 1, the routine
12 * returns the crudest estimate of the integral from A to B of f(x)
13 * dx. Subsequent calls with N=2,3,... (in that sequential order)
14 * will improve the accuracy by adding 2**(N-2) additional interior
17 * N.B., this function contains static state and IS NEITHER RENTRANT
22 trapzd (double (*func)(double),
26 long double x, tnm, sum, del;
32 it = 1; /* # of points to add on the next call */
33 s = 0.5 * (b - a) * (func(a) + func(b));
38 del = (b-a)/tnm; /* this is the spacing of the points to be added */
40 for (sum = 0.0, j = 1; j <= it; j++, x += del)
43 s = 0.5 * (s + (b-a) * sum/tnm); /* replace s by it's refined value */
49 * Returns the integral of the function FUNC from A to B. The
50 * parameters EPS can be set to the desired fractional accuracy and
51 * JMAX so that 2**(JMAX-1) is the maximum allowed number of steps.
52 * Integration is performed by Simpson's rule.
56 qsimp (double (*func)(double),
57 double a, /* lower limit */
58 double b) /* upper limit */
61 long double s, st, ost, os;
64 for (j = 1; j <= JMAX; j++){
65 st = trapzd (func, a, b, j);
66 s = (4.0 * st - ost)/3.0;
67 if (fabs (s - os) < EPS * fabs(os))
72 fprintf (stderr, "Too many steps in routine QSIMP\n");