! fvn comment : ! Modified version of the dqk15 quadpack routine from http://www.netlib.org/quadpack ! ! + The external 'f' function is a 2 parameters function f(x,y). The routine ! takes one more parameter 'x' and evaluate the integral of f against y between a and b ! for a given x subroutine dqk15_2d_inner(f,x,a,b,result,abserr,resabs,resasc) !***begin prologue dqk15 !***date written 800101 (yymmdd) !***revision date 830518 (yymmdd) !***category no. h2a1a2 !***keywords 15-point gauss-kronrod rules !***author piessens,robert,appl. math. & progr. div. - k.u.leuven ! de doncker,elise,appl. math. & progr. div - k.u.leuven !***purpose to compute i = integral of f over (a,b), with error ! estimate ! j = integral of abs(f) over (a,b) !***description ! ! integration rules ! standard fortran subroutine ! double precision version ! ! parameters ! on entry ! f - double precision ! function subprogram defining the integrand ! function f(x). the actual name for f needs to be ! declared e x t e r n a l in the calling program. ! ! a - double precision ! lower limit of integration ! ! b - double precision ! upper limit of integration ! ! on return ! result - double precision ! approximation to the integral i ! result is computed by applying the 15-point ! kronrod rule (resk) obtained by optimal addition ! of abscissae to the7-point gauss rule(resg). ! ! abserr - double precision ! estimate of the modulus of the absolute error, ! which should not exceed abs(i-result) ! ! resabs - double precision ! approximation to the integral j ! ! resasc - double precision ! approximation to the integral of abs(f-i/(b-a)) ! over (a,b) ! !***references (none) !***routines called d1mach !***end prologue dqk15 ! double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1, & epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc, & resg,resk,reskh,result,uflow,wg,wgk,xgk,x integer j,jtw,jtwm1 external f ! dimension fv1(7),fv2(7),wg(4),wgk(8),xgk(8) ! ! the abscissae and weights are given for the interval (-1,1). ! because of symmetry only the positive abscissae and their ! corresponding weights are given. ! ! xgk - abscissae of the 15-point kronrod rule ! xgk(2), xgk(4), ... abscissae of the 7-point ! gauss rule ! xgk(1), xgk(3), ... abscissae which are optimally ! added to the 7-point gauss rule ! ! wgk - weights of the 15-point kronrod rule ! ! wg - weights of the 7-point gauss rule ! ! ! gauss quadrature weights and kronron quadrature abscissae and weights ! as evaluated with 80 decimal digit arithmetic by l. w. fullerton, ! bell labs, nov. 1981. ! data wg ( 1) / 0.129484966168869693270611432679082d0 / data wg ( 2) / 0.279705391489276667901467771423780d0 / data wg ( 3) / 0.381830050505118944950369775488975d0 / data wg ( 4) / 0.417959183673469387755102040816327d0 / ! data xgk ( 1) / 0.991455371120812639206854697526329d0 / data xgk ( 2) / 0.949107912342758524526189684047851d0 / data xgk ( 3) / 0.864864423359769072789712788640926d0 / data xgk ( 4) / 0.741531185599394439863864773280788d0 / data xgk ( 5) / 0.586087235467691130294144838258730d0 / data xgk ( 6) / 0.405845151377397166906606412076961d0 / data xgk ( 7) / 0.207784955007898467600689403773245d0 / data xgk ( 8) / 0.000000000000000000000000000000000d0 / ! data wgk ( 1) / 0.022935322010529224963732008058970d0 / data wgk ( 2) / 0.063092092629978553290700663189204d0 / data wgk ( 3) / 0.104790010322250183839876322541518d0 / data wgk ( 4) / 0.140653259715525918745189590510238d0 / data wgk ( 5) / 0.169004726639267902826583426598550d0 / data wgk ( 6) / 0.190350578064785409913256402421014d0 / data wgk ( 7) / 0.204432940075298892414161999234649d0 / data wgk ( 8) / 0.209482141084727828012999174891714d0 / ! ! ! list of major variables ! ----------------------- ! ! centr - mid point of the interval ! hlgth - half-length of the interval ! absc - abscissa ! fval* - function value ! resg - result of the 7-point gauss formula ! resk - result of the 15-point kronrod formula ! reskh - approximation to the mean value of f over (a,b), ! i.e. to i/(b-a) ! ! machine dependent constants ! --------------------------- ! ! epmach is the largest relative spacing. ! uflow is the smallest positive magnitude. ! !***first executable statement dqk15 epmach = d1mach(4) uflow = d1mach(1) ! centr = 0.5d+00*(a+b) hlgth = 0.5d+00*(b-a) dhlgth = dabs(hlgth) ! ! compute the 15-point kronrod approximation to ! the integral, and estimate the absolute error. ! fc = f(x,centr) resg = fc*wg(4) resk = fc*wgk(8) resabs = dabs(resk) do 10 j=1,3 jtw = j*2 absc = hlgth*xgk(jtw) fval1 = f(x,centr-absc) fval2 = f(x,centr+absc) fv1(jtw) = fval1 fv2(jtw) = fval2 fsum = fval1+fval2 resg = resg+wg(j)*fsum resk = resk+wgk(jtw)*fsum resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2)) 10 continue do 15 j = 1,4 jtwm1 = j*2-1 absc = hlgth*xgk(jtwm1) fval1 = f(x,centr-absc) fval2 = f(x,centr+absc) fv1(jtwm1) = fval1 fv2(jtwm1) = fval2 fsum = fval1+fval2 resk = resk+wgk(jtwm1)*fsum resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2)) 15 continue reskh = resk*0.5d+00 resasc = wgk(8)*dabs(fc-reskh) do 20 j=1,7 resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh)) 20 continue result = resk*hlgth resabs = resabs*dhlgth resasc = resasc*dhlgth abserr = dabs((resk-resg)*hlgth) if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00) & abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00) if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1 & ((epmach*0.5d+02)*resabs,abserr) return end subroutine