k16pGaussLeg.f

Go to the documentation of this file.

!            test k16pGaussLeg.
!     implicit none
!     real*8 a, b, ans, func
!     integer n
!     external func
!     a = 0.
!     b = asin(1.d0) * 7.d0
!     do n = 2, 16
!        call k16pGaussLeg(func, a, b, n, ans)
!        write(*, *) ans
!     enddo
!     end
!     real*8 function func(x)
!     implicit none
!     real*8 x
!         func = sin(x)
!     end
!***************************************************
      subroutine k16pgaussleg(func, a, b, n, ans)
!***************************************************
!        Gauss-Legendre <= 16 points integration
!    quadrature. Adapted from Mori's code in
!    Maruzen book.
!
! func: real*8. input. integrand function name.
!    a: real*8. input. lower bound of  integration
!    b: real*8. input.  upper bound of integration 
!    n: integer. input  number of points
!                 to be used 2 <= n <= 16 
!  ans: real*8. output.  result of integration 
! -------------------------------------------------
      implicit none
      integer n
      real*8 func, a, b, ans
      external func

      real*8 c1, c2
      integer nh,i
      real*8 coef(0:16,16), weight(0:16,16)

      data (coef(i, 1),i=0, 0) /
     *   0.0000000000000000d+00 /
      data (weight(i, 1),i=0, 0) /
     *   0.2000000000000000d+01 /
      data (coef(i, 2),i=1, 1) /
     *   0.5773502691896257d+00 /
      data (weight(i, 2),i=1, 1) /
     *   0.9999999999999998d+00 /
      data (coef(i, 3),i=0, 1) /
     *   0.0000000000000000d+00,
     *   0.7745966692414832d+00 /
      data (weight(i, 3),i=0, 1) /
     *   0.8888888888888888d+00,
     *   0.5555555555555558d+00 /
      data (coef(i, 4),i=1, 2) /
     *   0.3399810435848563d+00,
     *   0.8611363115940525d+00 /
      data (weight(i, 4),i=1, 2) /
     *   0.6521451548625459d+00,
     *   0.3478548451374539d+00 /
      data (coef(i, 5),i=0, 2) /
     *   0.0000000000000000d+00,
     *   0.5384693101056831d+00,
     *   0.9061798459386639d+00 /
      data (weight(i, 5),i=0, 2) /
     *   0.5688888888888888d+00,
     *   0.4786286704993666d+00,
     *   0.2369268850561892d+00 /
      data (coef(i, 6),i=1, 3) /
     *   0.2386191860831969d+00,
     *   0.6612093864662644d+00,
     *   0.9324695142031520d+00 /
      data (weight(i, 6),i=1, 3) /
     *   0.4679139345726909d+00,
     *   0.3607615730481386d+00,
     *   0.1713244923791703d+00 /
      data (coef(i, 7),i=0, 3) /
     *   0.0000000000000000d+00,
     *   0.4058451513773972d+00,
     *   0.7415311855993944d+00,
     *   0.9491079123427584d+00 /
      data (weight(i, 7),i=0, 3) /
     *   0.4179591836734694d+00,
     *   0.3818300505051189d+00,
     *   0.2797053914892767d+00,
     *   0.1294849661688699d+00 /
      data (coef(i, 8),i=1, 4) /
     *   0.1834346424956498d+00,
     *   0.5255324099163289d+00,
     *   0.7966664774136267d+00,
     *   0.9602898564975361d+00 /
      data (weight(i, 8),i=1, 4) /
     *   0.3626837833783622d+00,
     *   0.3137066458778873d+00,
     *   0.2223810344533746d+00,
     *   0.1012285362903763d+00 /
      data (coef(i, 9),i=0, 4) /
     *   0.0000000000000000d+00,
     *   0.3242534234038089d+00,
     *   0.6133714327005903d+00,
     *   0.8360311073266357d+00,
     *   0.9681602395076260d+00 /
      data (weight(i, 9),i=0, 4) /
     *   0.3302393550012598d+00,
     *   0.3123470770400028d+00,
     *   0.2606106964029356d+00,
     *   0.1806481606948573d+00,
     *   0.8127438836157435d-01 /
      data (coef(i,10),i=1, 5) /
     *   0.1488743389816312d+00,
     *   0.4333953941292473d+00,
     *   0.6794095682990244d+00,
     *   0.8650633666889845d+00,
     *   0.9739065285171717d+00 /
      data (weight(i,10),i=1, 5) /
     *   0.2955242247147527d+00,
     *   0.2692667193099963d+00,
     *   0.2190863625159819d+00,
     *   0.1494513491505806d+00,
     *   0.6667134430868799d-01 /
      data (coef(i,11),i=0, 5) /
     *   0.0000000000000000d+00,
     *   0.2695431559523450d+00,
     *   0.5190961292068118d+00,
     *   0.7301520055740493d+00,
     *   0.8870625997680953d+00,
     *   0.9782286581460569d+00 /
      data (weight(i,11),i=0, 5) /
     *   0.2729250867779006d+00,
     *   0.2628045445102467d+00,
     *   0.2331937645919905d+00,
     *   0.1862902109277342d+00,
     *   0.1255803694649046d+00,
     *   0.5566856711617373d-01 /
      data (coef(i,12),i=1, 6) /
     *   0.1252334085114689d+00,
     *   0.3678314989981802d+00,
     *   0.5873179542866174d+00,
     *   0.7699026741943046d+00,
     *   0.9041172563704747d+00,
     *   0.9815606342467192d+00 /
      data (weight(i,12),i=1, 6) /
     *   0.2491470458134029d+00,
     *   0.2334925365383548d+00,
     *   0.2031674267230659d+00,
     *   0.1600783285433463d+00,
     *   0.1069393259953185d+00,
     *   0.4717533638651187d-01 /
      data (coef(i,13),i=0, 6) /
     *   0.0000000000000000d+00,
     *   0.2304583159551348d+00,
     *   0.4484927510364469d+00,
     *   0.6423493394403402d+00,
     *   0.8015780907333098d+00,
     *   0.9175983992229779d+00,
     *   0.9841830547185881d+00 /
      data (weight(i,13),i=0, 6) /
     *   0.2325515532308739d+00,
     *   0.2262831802628972d+00,
     *   0.2078160475368886d+00,
     *   0.1781459807619456d+00,
     *   0.1388735102197873d+00,
     *   0.9212149983772848d-01,
     *   0.4048400476531587d-01 /
      data (coef(i,14),i=1, 7) /
     *   0.1080549487073436d+00,
     *   0.3191123689278898d+00,
     *   0.5152486363581540d+00,
     *   0.6872929048116854d+00,
     *   0.8272013150697650d+00,
     *   0.9284348836635735d+00,
     *   0.9862838086968123d+00 /
      data (weight(i,14),i=1, 7) /
     *   0.2152638534631578d+00,
     *   0.2051984637212956d+00,
     *   0.1855383974779379d+00,
     *   0.1572031671581936d+00,
     *   0.1215185706879031d+00,
     *   0.8015808715976016d-01,
     *   0.3511946033175195d-01 /
      data (coef(i,15),i=0, 7) /
     *   0.0000000000000000d+00,
     *   0.2011940939974345d+00,
     *   0.3941513470775634d+00,
     *   0.5709721726085388d+00,
     *   0.7244177313601700d+00,
     *   0.8482065834104272d+00,
     *   0.9372733924007058d+00,
     *   0.9879925180204854d+00 /
      data (weight(i,15),i=0, 7) /
     *   0.2025782419255613d+00,
     *   0.1984314853271116d+00,
     *   0.1861610000155623d+00,
     *   0.1662692058169940d+00,
     *   0.1395706779261542d+00,
     *   0.1071592204671719d+00,
     *   0.7036604748810814d-01,
     *   0.3075324199611710d-01 /
      data (coef(i,16),i=1, 8) /
     *   0.9501250983763744d-01,
     *   0.2816035507792589d+00,
     *   0.4580167776572274d+00,
     *   0.6178762444026437d+00,
     *   0.7554044083550029d+00,
     *   0.8656312023878317d+00,
     *   0.9445750230732326d+00,
     *   0.9894009349916499d+00 /
      data (weight(i,16),i=1, 8) /
     *   0.1894506104550686d+00,
     *   0.1826034150449236d+00,
     *   0.1691565193950025d+00,
     *   0.1495959888165767d+00,
     *   0.1246289712555339d+00,
     *   0.9515851168249285d-01,
     *   0.6225352393864789d-01,
     *   0.2715245941175410d-01 /

      if( n .ge. 2 .and. n .le. 16) then
         c1 = (b + a) / 2
         c2 = (b - a) / 2

         if (mod(n,2) .eq. 0) then
            nh = n / 2
            ans = 0
         else
            nh = (n - 1) / 2
            ans = weight(0,n) * func(c1)
         endif

         do i = 1, nh
            ans = ans + 
     *           weight(i,n) * (func(c1 + c2 * coef(i,n))
     *          + func(c1 - c2 * coef(i,n)))
         enddo

         ans = c2 * ans
      else
         write (6,
     *  '(" k16pGaussLeg: n invalid=",i10)')
     *   n
         ans = 0.
      endif
      end