c---------------------------------------------------------------------- c This routine returns the result of numerically integrating a c polynomial of a given "size" over a given interval using Simpson's c method. c Arguments are as follows. c XLEFT and XRIGHT, the REAL endpoints of the interval; c NINTVL, the INTEGER number of subintervals; c F, the EXTERNAL REAL three-argument polynomial c evaluation function (F's arguments are the poly array, c its size, and the x value at which the polynomial is to c be evaluated); c POLY, the REAL array containing polynomial coefficients; c SIZE, the size of the array; c ERROR, a LOGICAL error indicator; and c AREA, a REAL variable in which PSIMP will return the desired c approximation to the value of the integral. subroutine psimp (xleft, xright, nintvl, . F, poly, size, error, area) real xleft, xright, area integer nintvl external F real poly (*), F integer size logical error real sumevn, sumodd, h, x integer k if (nintvl .lt. 1) then error = .true. return end if error = .false. h = (xright - xleft) / (2.0 * nintvl) k = 2 sumevn = 0.0 sumodd = F(poly, size, xleft+h) 1 if (k .le. 2 * nintvl - 2) then x = xleft + k*h sumevn = sumevn + F(poly, size, x) sumodd = sumodd + F(poly, size, x+h) k = k + 2 go to 1 endif area = (F(poly, size, xleft) + F(poly, size, xright) . + 4.0*sumodd + 2.0*sumevn) * h / 3.0 return end