double precision function rombint(f,a,b,tol,v1) c Rombint returns the integral from a to b of using Romberg integration. c The method converges provided that f(x) is continuous in (a,b). c f must be double precision and must be declared external in the calling c routine. tol indicates the desired relative accuracy in the integral. c c 1 free variable, v1, can be passed to the function. c parameter (MAXITER=40,MAXJ=5) implicit double precision (a-h,o-z) dimension g(MAXJ+1) double precision f external f c h=0.5d0*(b-a) gmax=h*(f(a,v1)+f(b,v1)) g(1)=gmax nint=1 error=1.0d20 i=0 10 i=i+1 if (i.gt.MAXITER.or.(i.gt.5.and.abs(error).lt.tol)) 2 go to 40 c Calculate next trapezoidal rule approximation to integral. g0=0.0d0 do 20 k=1,nint g0=g0+f(a+(k+k-1)*h,v1) 20 continue g0=0.5d0*g(1)+h*g0 h=0.5d0*h nint=nint+nint jmax=min(i,MAXJ) fourj=1.0d0 do 30 j=1,jmax c Use Richardson extrapolation. fourj=4.0d0*fourj g1=g0+(g0-g(j))/(fourj-1.0d0) g(j)=g0 g0=g1 30 continue if (abs(g0).gt.tol) then error=1.0d0-gmax/g0 else error=gmax end if gmax=g0 g(jmax+1)=g0 go to 10 40 rombint=g0 if (i.gt.MAXITER.and.abs(error).gt.tol) 2 write(*,*) 'Rombint failed to converge; integral, error=', 3 rombint,error return end