function arc_sine ( s ) !*****************************************************************************80 ! !! ARC_SINE computes the arc sine function, with argument truncation. ! ! Discussion: ! ! If you call your system ASIN routine with an input argument that is ! even slightly outside the range [-1.0, 1.0 ], you may get an unpleasant ! surprise (I did). ! ! This routine simply truncates arguments outside the range. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 28 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) S, the argument. ! ! Output, real ( kind = 8 ) ARC_SINE, an angle whose sine is S. ! implicit none real ( kind = 8 ) arc_sine real ( kind = 8 ) s real ( kind = 8 ) s2 s2 = s s2 = max ( s2, -1.0D+00 ) s2 = min ( s2, +1.0D+00 ) arc_sine = asin ( s2 ) return end subroutine ball_f1_nd ( func, n, center, r, result ) !*****************************************************************************80 ! !! BALL_F1_ND approximates an integral inside a ball in ND. ! ! Integration region: ! ! sum ( X(1:N) - CENTER(1:N) )^2 <= R * R. ! ! Discussion: ! ! An (N+1)*2**N point 5-th degree formula is used, Stroud number SN:5-6. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 December 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates F at the N-vector X, of the form ! function func ( n, x ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) func ! real ( kind = 8 ) x(n) ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Input, real ( kind = 8 ) CENTER(N), the center of the ball. ! ! Input, real ( kind = 8 ) R, the radius of the ball. ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) ball_volume_nd real ( kind = 8 ) center(n) real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) ihi integer ( kind = 4 ) itemp integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) khi integer ( kind = 4 ) ktemp real ( kind = 8 ) quad real ( kind = 8 ) r real ( kind = 8 ) result real ( kind = 8 ) t real ( kind = 8 ) temp real ( kind = 8 ) u real ( kind = 8 ) u2 real ( kind = 8 ) v real ( kind = 8 ) volume real ( kind = 8 ) w real ( kind = 8 ) x(n) real ( kind = 8 ) y if ( r == 0.0D+00 ) then result = 0.0D+00 return end if u2 = ( 1.0D+00 - 2.0D+00 * sqrt ( 1.0D+00 / real ( n + 4, kind = 8 ) ) ) & / real ( n + 2, kind = 8 ) u = sqrt ( u2 ) x(1:n) = center(1:n) - r * u w = 1.0D+00 / real ( ( n + 1 ) * 2**n, kind = 8 ) quad = 0.0D+00 ihi = 2**n do i = 1, ihi itemp = i - 1 do j = 1, n u = ( center(j) - x(j) ) / r if ( mod ( itemp, 2 ) == 1 ) then x(j) = center(j) - abs ( x(j) - center(j) ) else x(j) = center(j) + abs ( x(j) - center(j) ) end if itemp = itemp / 2 end do quad = quad + w * func ( n, x ) end do temp = sqrt ( real ( n + 4, kind = 8 ) ) t = sqrt ( 2.0D+00 * real ( n + 1, kind = 8 ) / real ( n + 2, kind = 8 ) ) & / ( real ( n, kind = 8 ) * temp ) y = ( 1.0D+00 + 2.0D+00 / ( real ( n, kind = 8 ) * temp ) ) & / real ( n + 2, kind = 8 ) v = sqrt ( y - t ) u = sqrt ( y + real ( n - 1, kind = 8 ) * t ) khi = 2**n do i = 1, n x(1:n) = center(1:n) - r * v x(i) = center(i) - r * u do k = 1, khi ktemp = k - 1 do j = 1, n if ( mod ( ktemp, 2 ) == 1 ) then x(j) = center(j) - abs ( x(j) - center(j) ) else x(j) = center(j) + abs ( x(j) - center(j) ) end if ktemp = ktemp / 2 end do quad = quad + w * func ( n, x ) end do x(i) = center(i) - r * v end do volume = ball_volume_nd ( n, r ) result = quad * volume return end subroutine ball_f3_nd ( func, n, center, r, result ) !*****************************************************************************80 ! !! BALL_F3_ND approximates an integral inside a ball in ND. ! ! Integration region: ! ! sum ( X(1:N) - CENTER(1:N) )^2 <= R * R. ! ! Discussion: ! ! A 2**(N+1)-1 point 5-th degree formula is used, Stroud number SN:5-4. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 21 November 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates F at the N-vector X, of the form ! function func ( n, x ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) func ! real ( kind = 8 ) x(n) ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Input, real ( kind = 8 ) CENTER(N), the center of the ball. ! ! Input, real ( kind = 8 ) R, the radius of the ball. ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) ball_volume_nd real ( kind = 8 ) center(n) real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) jtemp integer ( kind = 4 ) k real ( kind = 8 ) quad real ( kind = 8 ) r real ( kind = 8 ) result real ( kind = 8 ) ri real ( kind = 8 ) s real ( kind = 8 ) volume real ( kind = 8 ) weight real ( kind = 8 ) x(n) if ( r == 0.0D+00 ) then result = 0.0D+00 return end if quad = 0.0D+00 ! ! The first point is the center of the ball. ! x(1:n) = center(1:n) weight = 4.0D+00 / real ( ( n + 2 ) * ( n + 2 ), kind = 8 ) quad = quad + weight * func ( n, x ) s = 1.0D+00 / sqrt ( real ( n + 4, kind = 8 ) ) do i = 1, n ri = sqrt ( real ( i + 2, kind = 8 ) / real ( n + 4, kind = 8 ) ) ! ! Set up the first point, with (I-1) zeroes, RI, and then N-I S's. ! do j = 1, n if ( j < i ) then x(j) = center(j) else if ( j == i ) then x(j) = center(j) + r * ri else x(j) = center(j) + r * s end if end do weight = 2.0D+00**( i - n ) * real ( n + 4, kind = 8 ) & / real ( ( i + 1 ) * ( i + 2 ) * ( n + 2 ), kind = 8 ) ! ! Now go through all sign permutations of the basic point. ! do j = 1, 2**(n+1-i) jtemp = j - 1 do k = i, n if ( mod ( jtemp, 2 ) == 1 ) then x(k) = center(k) - abs ( x(k) - center(k) ) else x(k) = center(k) + abs ( x(k) - center(k) ) end if jtemp = jtemp / 2 end do quad = quad + weight * func ( n, x ) end do end do volume = ball_volume_nd ( n, r ) result = quad * volume return end function ball_monomial_nd ( n, p, r ) !*****************************************************************************80 ! !! BALL_MONOMIAL_ND integrates a monomial on a ball in ND. ! ! Integration region: ! ! sum ( X(1:N)^2 ) <= R * R ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 31 May 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Gerald Folland, ! How to Integrate a Polynomial Over a Sphere, ! American Mathematical Monthly, ! Volume 108, Number 5, May 2001, pages 446-448. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Input, integer ( kind = 4 ) P(N), the exponents of X(1) through X(N) in the monomial. ! The exponents P(N) must be nonnegative. ! ! Input, real ( kind = 8 ) R, the radius of the ball. ! ! Output, real ( kind = 8 ) BALL_MONOMIAL_ND, the integral of ! X1**P(1)*X2**P(2)*...*XN**P(N) over the ball. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) ball_monomial_nd integer ( kind = 4 ) p(n) real ( kind = 8 ) power real ( kind = 8 ) r real ( kind = 8 ) sphere_unit_monomial_nd power = real ( sum ( p ) + n, kind = 8 ) ball_monomial_nd = sphere_unit_monomial_nd ( n, p ) * r**power / power return end subroutine ball_unit_07_3d ( func, result ) !*****************************************************************************80 ! !! BALL_UNIT_07_3D approximates an integral inside the unit ball in 3D. ! ! Integration region: ! ! X*X + Y*Y + Z*Z <= 1. ! ! Discussion: ! ! A 64 point 7-th degree formula is used. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 21 November 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates F(X,Y,Z), of the form ! function func ( x, y, z ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! real ( kind = 8 ) z ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ), parameter :: order = 4 real ( kind = 8 ) angle real ( kind = 8 ) ball_unit_volume_3d real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) k real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) quad real ( kind = 8 ) result real ( kind = 8 ) volume real ( kind = 8 ) w real ( kind = 8 ) weight1(order) real ( kind = 8 ) weight2(order) real ( kind = 8 ) weight3(order) real ( kind = 8 ) x real ( kind = 8 ) xtab1(order) real ( kind = 8 ) xtab2(order) real ( kind = 8 ) xtab3(order) real ( kind = 8 ) y real ( kind = 8 ) z ! ! This is the 5 point Gauss-Legendre rule, ! but with the midpoint deleted, and with different weights. ! xtab1(1:4) = (/ & -0.906179845938663992797626878299D+00, & -0.538469310105683091036314420700D+00, & 0.538469310105683091036314420700D+00, & 0.906179845938663992797626878299D+00 /) weight1(1:4) = (/ & 0.19455533421780251826D+00, & 0.13877799911553081506D+00, & 0.13877799911553081506D+00, & 0.19455533421780251826D+00 /) ! ! Set XTAB2 and WEIGHT2. ! do j = 1, order angle = pi * real ( 2 * j - 1, kind = 8 ) & / real ( 2 * order, kind = 8 ) xtab2(j) = cos ( angle ) end do weight2(1:order) = 1.0D+00 ! ! Set XTAB3 and WEIGHT3 for the interval [-1,1]. ! call legendre_set ( order, xtab3, weight3 ) w = 3.0D+00 / 16.0D+00 quad = 0.0D+00 do i = 1, order do j = 1, order do k = 1, order x = xtab1(i) * sqrt ( 1.0D+00 - xtab2(j) * xtab2(j) ) & * sqrt ( 1.0D+00 - xtab3(k) * xtab3(k) ) y = xtab1(i) * xtab2(j) * sqrt ( 1.0D+00 - xtab3(k) * xtab3(k) ) z = xtab1(i) * xtab3(k) quad = quad + w * weight1(i) * weight2(j) * weight3(k) & * func ( x, y, z ) end do end do end do volume = ball_unit_volume_3d ( ) result = quad * volume return end subroutine ball_unit_14_3d ( func, result ) !*****************************************************************************80 ! !! BALL_UNIT_14_3D approximates an integral inside the unit ball in 3D. ! ! Integration region: ! ! X*X + Y*Y + Z*Z <= 1. ! ! Discussion: ! ! A 288 point 14-th degree formula is used, Stroud number S3:14-1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 21 November 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates F(X,Y,Z), of the form ! function func ( x, y, z ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! real ( kind = 8 ) z ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none real ( kind = 8 ) ball_unit_volume_3d real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) quad real ( kind = 8 ), save, dimension ( 4 ) :: r = (/ & 0.968160240D+00, 0.836031107D+00, 0.613371433D+00, 0.324253423D+00 /) real ( kind = 8 ) result real ( kind = 8 ) temp real ( kind = 8 ) volume real ( kind = 8 ) w1 real ( kind = 8 ) w2 real ( kind = 8 ), save, dimension ( 4 ) :: weight = (/ & 0.076181268D+00, 0.126263673D+00, 0.098048133D+00, 0.032840260D+00 /) real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( 5 ) :: xtab = (/ & -0.151108275D+00, 0.315838353D+00, 0.346307112D+00, & -0.101808787D+00, -0.409228403D+00 /) real ( kind = 8 ) y real ( kind = 8 ), save, dimension ( 5 ) :: ytab = (/ & 0.155240600D+00, 0.257049387D+00, 0.666277790D+00, & 0.817386065D+00, 0.501547712D+00 /) real ( kind = 8 ) z real ( kind = 8 ), save, dimension ( 5 ) :: ztab = (/ & 0.976251323D+00, 0.913330032D+00, 0.660412970D+00, & 0.567022920D+00, 0.762221757D+00 /) quad = 0.0D+00 do m = 1, 4 w1 = 125.0D+00 * weight(m) / 3360.0D+00 x = 0.525731112D+00 * r(m) y = 0.850650808D+00 * r(m) z = 0.0D+00 do j = 1, 2 x = -x do k = 1, 2 y = -y do l = 1, 3 call r8_swap3 ( x, y, z ) quad = quad + w1 * func ( x, y, z ) end do end do end do w2 = 143.0D+00 * weight(m) / 3360.0D+00 do n = 1, 5 x = xtab(n) * r(m) y = ytab(n) * r(m) z = ztab(n) * r(m) do i = 1, 3 temp = x x = z z = -y y = -temp do j = 1, 3 call r8_swap3 ( x, y, z ) quad = quad + w2 * func ( x, y, z ) end do y = -y z = -z quad = quad + w2 * func ( x, y, z ) end do end do end do volume = ball_unit_volume_3d ( ) result = quad * volume return end subroutine ball_unit_15_3d ( func, result ) !*****************************************************************************80 ! !! BALL_UNIT_15_3D approximates an integral inside the unit ball in 3D. ! ! Integration region: ! ! X*X + Y*Y + Z*Z <= 1. ! ! Discussion: ! ! A 512 point 15-th degree formula is used. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 28 October 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates F(X,Y,Z), of the form ! function func ( x, y, z ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! real ( kind = 8 ) z ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ), parameter :: order1 = 4 integer ( kind = 4 ), parameter :: order2 = 8 real ( kind = 8 ) ball_unit_volume_3d real ( kind = 8 ) cj real ( kind = 8 ) ck real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) k real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) quad real ( kind = 8 ) result real ( kind = 8 ) sj real ( kind = 8 ) sk real ( kind = 8 ) volume real ( kind = 8 ) w real ( kind = 8 ), save, dimension ( order1 ) :: weight1 = (/ & 0.0328402599D+00, 0.0980481327D+00, 0.1262636728D+00, 0.0761812678D+00 /) real ( kind = 8 ) weight2(order2) real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( order1 ) :: xtab1 = (/ & 0.3242534234D+00, 0.6133714327D+00, 0.8360311073D+00, 0.9681602395D+00 /) real ( kind = 8 ) xtab2(order2) real ( kind = 8 ) y real ( kind = 8 ) z call legendre_set ( order2, xtab2, weight2 ) w = 3.0D+00 / 32.0D+00 quad = 0.0D+00 do i = 1, order1 do j = 1, order2 sj = xtab2(j) cj = sqrt ( 1.0D+00 - sj * sj ) do k = 1, 16 sk = sin ( real ( k, kind = 8 ) * pi / 8.0D+00 ) ck = cos ( real ( k, kind = 8 ) * pi / 8.0D+00 ) x = xtab1(i) * cj * ck y = xtab1(i) * cj * sk z = xtab1(i) * sj quad = quad + w * weight1(i) * weight2(j) * func ( x, y, z ) end do end do end do volume = ball_unit_volume_3d ( ) result = quad * volume return end subroutine ball_unit_f1_nd ( func, n, result ) !*****************************************************************************80 ! !! BALL_UNIT_F1_ND approximates an integral inside the unit ball in ND. ! ! Integration region: ! ! sum ( X(1:N)^2 ) <= 1. ! ! Discussion: ! ! An (N+1)*2**N point 5-th degree formula is used, Stroud number SN:5-6. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 21 November 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates F at the N-vector X, of the form ! function func ( n, x ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) func ! real ( kind = 8 ) x(n) ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) ball_unit_volume_nd real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) ihi integer ( kind = 4 ) itemp integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) khi integer ( kind = 4 ) ktemp real ( kind = 8 ) quad real ( kind = 8 ) result real ( kind = 8 ) t real ( kind = 8 ) temp real ( kind = 8 ) u real ( kind = 8 ) u2 real ( kind = 8 ) v real ( kind = 8 ) volume real ( kind = 8 ) w real ( kind = 8 ) x(n) real ( kind = 8 ) y u2 = ( 1.0D+00 - 2.0D+00 * sqrt ( 1.0D+00 / real ( n + 4, kind = 8 ) ) ) & / real ( n + 2, kind = 8 ) u = sqrt ( u2 ) x(1:n) = -u w = 1.0D+00 / real ( ( n + 1 ) * 2**n, kind = 8 ) quad = 0.0D+00 ihi = 2**n do i = 1, ihi itemp = i - 1 do j = 1, n if ( mod ( itemp, 2 ) == 1 ) then x(j) = -abs ( x(j) ) else x(j) = abs ( x(j) ) end if itemp = itemp / 2 end do quad = quad + w * func ( n, x ) end do temp = sqrt ( real ( n + 4, kind = 8 ) ) t = sqrt ( 2.0D+00 * real ( n + 1, kind = 8 ) / real ( n + 2, kind = 8 ) ) & / ( real ( n, kind = 8 ) * temp ) y = ( 1.0D+00 + 2.0D+00 / ( real ( n, kind = 8 ) * temp ) ) & / real ( n + 2, kind = 8 ) v = sqrt ( y - t ) u = sqrt ( y + real ( n - 1, kind = 8 ) * t ) khi = 2**n do i = 1, n x(1:n) = -v x(i) = -u do k = 1, khi ktemp = k - 1 do j = 1, n if ( mod ( ktemp, 2 ) == 1 ) then x(j) = -abs ( x(j) ) else x(j) = abs ( x(j) ) end if ktemp = ktemp / 2 end do quad = quad + w * func ( n, x ) end do x(i) = -v end do volume = ball_unit_volume_nd ( n ) result = quad * volume return end subroutine ball_unit_f3_nd ( func, n, result ) !*****************************************************************************80 ! !! BALL_UNIT_F3_ND approximates an integral inside the unit ball in ND. ! ! Integration region: ! ! sum ( X(1:N)^2 ) <= 1. ! ! Discussion: ! ! A 2**(N+1)-1 point 5-th degree formula is used, Stroud number SN:5-4. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 21 November 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates F at the N-vector X, of the form ! function func ( n, x ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) func ! real ( kind = 8 ) x(n) ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) ball_unit_volume_nd real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) jtemp integer ( kind = 4 ) k real ( kind = 8 ) quad real ( kind = 8 ) result real ( kind = 8 ) ri real ( kind = 8 ) s real ( kind = 8 ) volume real ( kind = 8 ) weight real ( kind = 8 ) x(n) quad = 0.0D+00 ! ! The first point is the center of the ball. ! x(1:n) = 0.0D+00 weight = 4.0D+00 / real ( ( n + 2 ) * ( n + 2 ), kind = 8 ) quad = quad + weight * func ( n, x ) s = 1.0D+00 / sqrt ( real ( n + 4, kind = 8 ) ) do i = 1, n ri = sqrt ( real ( i + 2, kind = 8 ) / real ( n + 4, kind = 8 ) ) ! ! Set up the first point, with (I-1) zeroes, RI, and then N-I S's. ! do j = 1, n if ( j < i ) then x(j) = 0.0D+00 else if ( j == i ) then x(j) = ri else x(j) = s end if end do weight = 2.0D+00**( i - n ) * real ( n + 4, kind = 8 ) & / real ( ( i + 1 ) * ( i + 2 ) * ( n + 2 ), kind = 8 ) ! ! Now go through all sign permutations of the basic point. ! do j = 1, 2**(n+1-i) jtemp = j - 1 do k = i, n if ( mod ( jtemp, 2 ) == 1 ) then x(k) = -abs ( x(k) ) else x(k) = abs ( x(k) ) end if jtemp = jtemp / 2 end do quad = quad + weight * func ( n, x ) end do end do volume = ball_unit_volume_nd ( n ) result = quad * volume return end function ball_unit_volume_3d ( ) !*****************************************************************************80 ! !! BALL_UNIT_VOLUME_3D computes the volume of the unit ball in 3D. ! ! Integration region: ! ! X*X + Y*Y + Z*Z <= 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 20 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) BALL_UNIT_VOLUME_3D, the volume of the ball. ! implicit none real ( kind = 8 ) ball_unit_volume_3d real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 ball_unit_volume_3d = ( 4.0D+00 / 3.0D+00 ) * pi return end function ball_unit_volume_nd ( n ) !*****************************************************************************80 ! !! BALL_UNIT_VOLUME_ND computes the volume of the unit ball in ND. ! ! Integration region: ! ! sum ( X(1:N)^2 ) <= 1. ! ! Discussion: ! ! N Volume ! ! 2 PI ! 3 (4/3) * PI ! 4 (1/2) * PI^2 ! 5 (8/15) * PI^2 ! 6 (1/6) * PI^3 ! 7 (16/105) * PI^3 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Output, real ( kind = 8 ) BALL_UNIT_VOLUME_ND, the volume of the ball. ! implicit none real ( kind = 8 ) ball_unit_volume_nd integer ( kind = 4 ) i integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) volume if ( mod ( n, 2 ) == 0 ) then m = n / 2 volume = ( pi )**m do i = 1, m volume = volume / real ( i, kind = 8 ) end do else m = ( n - 1 ) / 2 volume = ( pi )**m * 2.0D+00**n do i = m+1, 2*m+1 volume = volume / real ( i, kind = 8 ) end do end if ball_unit_volume_nd = volume return end function ball_volume_3d ( r ) !*****************************************************************************80 ! !! BALL_VOLUME_3D computes the volume of a ball in 3D. ! ! Integration region: ! ! X*X + Y*Y + Z*Z <= R * R ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R, the radius of the ball. ! ! Output, real ( kind = 8 ) BALL_VOLUME_3D, the volume of the ball. ! implicit none real ( kind = 8 ) ball_volume_3d real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) r ball_volume_3d = ( 4.0D+00 / 3.0D+00 ) * pi * r**3 return end function ball_volume_nd ( n, r ) !*****************************************************************************80 ! !! BALL_VOLUME_ND computes the volume of a ball in ND. ! ! Integration region: ! ! sum ( X(1:N)^2 ) <= R * R ! ! Discussion: ! ! N Volume ! ! 2 PI * R^2 ! 3 (4/3) * PI * R^3 ! 4 (1/2) * PI^2 * R^4 ! 5 (8/15) * PI^2 * R^5 ! 6 (1/6) * PI^3 * R^6 ! 7 (16/105) * PI^3 * R^7 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Input, real ( kind = 8 ) R, the radius of the ball. ! ! Output, real ( kind = 8 ) BALL_VOLUME_ND, the volume of the ball. ! implicit none real ( kind = 8 ) ball_unit_volume_nd real ( kind = 8 ) ball_volume_nd integer ( kind = 4 ) n real ( kind = 8 ) r ball_volume_nd = ball_unit_volume_nd ( n ) * r**n return end subroutine circle_annulus ( func, center, radius1, radius2, nr, result ) !*****************************************************************************80 ! !! CIRCLE_ANNULUS approximates an integral in an annulus. ! ! Integration region: ! ! RADIUS1^2 <= ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 <= RADIUS2^2 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 January 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! William Peirce, ! Numerical Integration Over the Planar Annulus, ! Journal of the Society for Industrial and Applied Mathematics, ! Volume 5, Number 2, June 1957, pages 66-73. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied function of two ! variables which is to be integrated, of the form: ! function func ( x, y ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! ! Input, real ( kind = 8 ) CENTER(2), the center of the circle. ! ! Input, real ( kind = 8 ) RADIUS1, RADIUS2, the radii of the circles. ! ! Input, integer ( kind = 4 ) NR, the order of the rule. This quantity specifies ! the number of distinct radii to use. The number of angles used will ! be 4*NR, for a total of 4*NR*NR points. ! ! Output, real ( kind = 8 ) RESULT, the approximation to the integral. ! implicit none integer ( kind = 4 ) nr integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) a real ( kind = 8 ) area real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) center(dim_num) real ( kind = 8 ) circle_annulus_area_2d real ( kind = 8 ) d real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) nt real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) quad real ( kind = 8 ) ra(nr) real ( kind = 8 ) radius1 real ( kind = 8 ) radius2 real ( kind = 8 ) result real ( kind = 8 ) rw(nr) real ( kind = 8 ) t real ( kind = 8 ) tw real ( kind = 8 ) x real ( kind = 8 ) y ! ! Choose radial abscissas and weights. ! call legendre_set ( nr, ra, rw ) a = -1.0D+00 b = +1.0D+00 c = radius1 * radius1 d = radius2 * radius2 call rule_adjust ( a, b, c, d, nr, ra, rw ) ra(1:nr) = sqrt ( ra(1:nr) ) rw(1:nr) = rw(1:nr) / ( radius2 - radius1 ) / ( radius2 + radius1 ) ! ! Set angular abscissas and weights. ! nt = 4 * nr tw = 1.0D+00 / real ( nt, kind = 8 ) ! ! Approximate the integral. ! quad = 0.0D+00 do i = 1, nt t = 2.0D+00 * pi * real ( i - 1, kind = 8 ) / real ( nt, kind = 8 ) do j = 1, nr x = center(1) + ra(j) * cos ( t ) y = center(2) + ra(j) * sin ( t ) quad = quad + tw * rw(j) * func ( x, y ) end do end do area = circle_annulus_area_2d ( radius1, radius2 ) result = quad * area return end function circle_annulus_area_2d ( radius1, radius2 ) !*****************************************************************************80 ! !! CIRCLE_ANNULUS_AREA_2D returns the area of a circular annulus in 2D. ! ! Integration region: ! ! RADIUS1^2 <= ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 <= RADIUS2^2 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) RADIUS1, RADIUS2, the radii of the circles. ! ! Output, real ( kind = 8 ) CIRCLE_ANNULUS_AREA_2D, the area of the annulus. ! implicit none real ( kind = 8 ) circle_annulus_area_2d real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) radius1 real ( kind = 8 ) radius2 circle_annulus_area_2d = pi * ( radius1 + radius2 ) & * ( radius2 - radius1 ) return end subroutine circle_annulus_sector ( func, center, radius1, radius2, theta1, & theta2, nr, result ) !*****************************************************************************80 ! !! CIRCLE_ANNULUS_SECTOR approximates an integral in a circular annulus sector. ! ! Discussion: ! ! A circular annulus sector comprises the area between two concentric ! circles and two concentric rays. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 January 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! William Peirce, ! Numerical Integration Over the Planar Annulus, ! Journal of the Society for Industrial and Applied Mathematics, ! Volume 5, Number 2, June 1957, pages 66-73. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied function of two ! variables which is to be integrated, of the form: ! function func ( x, y ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! ! Input, real ( kind = 8 ) CENTER(2), the center of the circle. ! ! Input, real ( kind = 8 ) RADIUS1, RADIUS2, the radii of the circles. ! ! Input, real ( kind = 8 ) THETA1, THETA2, the angles defining the sector. ! The sector is measured from THETA1 to THETA2. ! ! Input, integer ( kind = 4 ) NR, the order of the rule. This quantity specifies ! the number of distinct radii to use. The number of angles used will ! be 4*NR, for a total of 4*NR*NR points. ! ! Output, real ( kind = 8 ) RESULT, the approximation to the integral. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 integer ( kind = 4 ) nr real ( kind = 8 ) a real ( kind = 8 ) area real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) center(dim_num) real ( kind = 8 ) circle_annulus_sector_area_2d real ( kind = 8 ) d real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) nt real ( kind = 8 ) quad real ( kind = 8 ) ra(nr) real ( kind = 8 ) radius1 real ( kind = 8 ) radius2 real ( kind = 8 ) result real ( kind = 8 ) rw(nr) real ( kind = 8 ) ta(4*nr) real ( kind = 8 ) theta1 real ( kind = 8 ) theta2 real ( kind = 8 ) tw(4*nr) real ( kind = 8 ) x real ( kind = 8 ) y ! ! Set the radial abscissas and weights. ! call legendre_set ( nr, ra, rw ) a = -1.0D+00 b = +1.0D+00 c = radius1 * radius1 d = radius2 * radius2 call rule_adjust ( a, b, c, d, nr, ra, rw ) ra(1:nr) = sqrt ( ra(1:nr) ) rw(1:nr) = rw(1:nr) / ( radius2 - radius1 ) / ( radius2 + radius1 ) ! ! Pick angles evenly spaced between THETA1 and THETA2, but do not ! include the endpoints, and use a half interval for the first and last. ! nt = 4 * nr call tvec_even_bracket3 ( nt, theta1, theta2, ta ) tw(1:nt) = 1.0D+00 / real ( nt, kind = 8 ) ! ! Approximate the integral. ! quad = 0.0D+00 do i = 1, nt do j = 1, nr x = center(1) + ra(j) * cos ( ta(i) ) y = center(2) + ra(j) * sin ( ta(i) ) quad = quad + tw(i) * rw(j) * func ( x, y ) end do end do area = circle_annulus_sector_area_2d ( radius1, radius2, theta1, theta2 ) result = quad * area return end function circle_annulus_sector_area_2d ( radius1, radius2, theta1, theta2 ) !*****************************************************************************80 ! !! CIRCLE_ANNULUS_SECTOR_AREA_2D returns the area of a circular annulus sector in 2D. ! ! Discussion: ! ! A circular annulus sector comprises the area between two concentric ! circles and two concentric rays. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) RADIUS1, RADIUS2, the radii of the circles. ! ! Input, real ( kind = 8 ) THETA1, THETA2, the angles of the rays. ! Ordinarily, (THETA2-THETA1) is between 0 and 2*PI. ! ! Output, real ( kind = 8 ) CIRCLE_ANNULUS_SECTOR_AREA_2D, the area of the ! circular annulus sector. ! implicit none real ( kind = 8 ) circle_annulus_sector_area_2d real ( kind = 8 ) radius1 real ( kind = 8 ) radius2 real ( kind = 8 ) theta1 real ( kind = 8 ) theta2 circle_annulus_sector_area_2d = 0.5D+00 * ( radius1 + radius2 ) & * ( radius2 - radius1 ) * ( theta2 - theta1 ) return end function circle_area_2d ( r ) !*****************************************************************************80 ! !! CIRCLE_AREA_2D returns the area of a circle in 2D. ! ! Integration region: ! ! ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 <= R * R ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R, the radius of the circle. ! ! Output, real ( kind = 8 ) CIRCLE_AREA_2D, the area of the circle. ! implicit none real ( kind = 8 ) circle_area_2d real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) r circle_area_2d = pi * r * r return end subroutine circle_cap_area_2d ( r, h, area ) !*****************************************************************************80 ! !! CIRCLE_CAP_AREA_2D computes the area of a circle cap in 2D. ! ! Discussion: ! ! Draw any radius R of the circle and denote as P the point where the ! radius intersects the circle. Now consider the point Q which lies ! on the radius and which is H units from P. The line which is ! perpendicular to the radius R and passes through Q divides the ! circle into two pieces. The piece including the point P is the ! circular cap of height (or thickness) H. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 May 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R, the radius of the circle. ! ! Input, real ( kind = 8 ) H, the "height" of the circle cap. ! ! Output, real ( kind = 8 ) AREA, the area of the circle cap. ! implicit none real ( kind = 8 ) arc_sine real ( kind = 8 ) area real ( kind = 8 ) h real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) r real ( kind = 8 ) theta if ( h <= 0.0D+00 ) then area = 0.0D+00 else if ( h <= r ) then theta = 2.0D+00 * arc_sine ( sqrt ( h * ( 2.0D+00 * r - h ) ) / r ) area = r * r * ( theta - sin ( theta ) ) / 2.0D+00 else if ( h <= 2.0D+00 * r ) then theta = 2.0D+00 * arc_sine ( sqrt ( h * ( 2.0D+00 * r - h ) ) / r ) area = r * r * ( pi - ( theta - sin ( theta ) ) / 2.0D+00 ) else if ( 2.0D+00 * r <= h ) then area = pi * r * r end if return end subroutine circle_cum ( func, center, radius, order, result ) !*****************************************************************************80 ! !! CIRCLE_CUM approximates an integral on the circumference of a circle in 2D. ! ! Integration region: ! ! ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 <= R * R ! ! Discussion: ! ! An ORDER point, (ORDER-1)-th degree formula is used, ! Stroud number U2:M-1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 September 1998 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied function of two ! variables which is to be integrated, of the form: ! function func ( x, y ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! ! Input, real ( kind = 8 ) CENTER(2), the coordinates of the center of ! the circle. ! ! Input, real ( kind = 8 ) RADIUS, the radius of the circle. ! ! Input, integer ( kind = 4 ) ORDER, the number of points to use. ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 real ( kind = 8 ) angle real ( kind = 8 ) center(dim_num) real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) order real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) quad real ( kind = 8 ) radius real ( kind = 8 ) result real ( kind = 8 ) volume real ( kind = 8 ) x real ( kind = 8 ) y quad = 0.0D+00 do i = 1, order angle = real ( 2 * i, kind = 8 ) * pi / real ( order, kind = 8 ) x = center(1) + radius * cos ( angle ) y = center(2) + radius * sin ( angle ) quad = quad + func ( x, y ) end do quad = quad / real ( order, kind = 8 ) volume = pi * radius * radius result = quad * volume return end subroutine circle_rt_set ( rule, nr, nt, nc, ra, rw, ta, tw, cw ) !*****************************************************************************80 ! !! CIRCLE_RT_SET sets an R, THETA product quadrature rule in the unit circle. ! ! Discussion: ! ! For a given value of RULE, here are the number of points used at the ! center (NC), the number of points along the radial direction (NR) and ! the number of points along the theta direction (NT). The total number ! of points in the rule will be ! ! Total = NC + NR * NT. ! ! The user, when choosing RULE, must allocate enough space in the arrays ! RA, RW, TA and TW for the resulting values of NR and NT. ! ! RULE NC NR NT Total ! ---- -- -- -- ----- ! 1 1 0 0 1 ! 2 0 1 4 4 ! 3 1 1 4 5 ! 4 1 1 6 7 ! 5 1 2 4 9 ! 6 0 3 4 12 ! 7 1 2 10 21 ! 8 0 4 16 64 ! 9 0 5 20 120 ! ! The integral of F(X,Y) over the unit circle is approximated by ! ! Integral ( X*X + Y*Y <= 1 ) F(X,Y) dx dy ! = Integral ( 0 <= R <= 1, 0 <= T <= 2PI ) F(R*cos(T),R*sin(T)) r dr dt ! = approximately ! CW * F(0,0) ! + sum ( 1 <= I <= NR ) Sum ( 1 <= J <= NT ) ! RW(I) * TW(J) * F ( R(I) * cos ( TA(J) ), R(I) * sin ( TA(J) ) ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 March 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964, ! ISBN: 0-486-61272-4, ! LC: QA47.A34. ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, integer ( kind = 4 ) RULE, the rule desired. ! ! Input, integer ( kind = 4 ) NR, the number of R abscissas. ! ! Input, integer ( kind = 4 ) NT, the number of Theta abscissas. ! ! Input, integer ( kind = 4 ) NC, the number of center abscissas (0 or 1 ). ! ! Output, real ( kind = 8 ) RA(NR), RW(NR), the R abscissas and weights. ! ! Output, real ( kind = 8 ) TA(NT), TW(NT), the THETA abscissas and weights. ! ! Output, real ( kind = 8 ) ZW, the weight to use for the center. ! implicit none integer ( kind = 4 ) nr integer ( kind = 4 ) nt real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) cw real ( kind = 8 ) d integer ( kind = 4 ) nc real ( kind = 8 ) ra(nr) real ( kind = 8 ) rw(nr) integer ( kind = 4 ) rule real ( kind = 8 ) ta(nt) real ( kind = 8 ) tw(nt) real ( kind = 8 ) u real ( kind = 8 ) v real ( kind = 8 ) w if ( rule == 1 ) then cw = 1.0D+00 else if ( rule == 2 ) then ra(1) = 0.5D+00 rw(1) = 1.0D+00 call tvec_even2 ( nt, ta ) tw(1:nt) = 1.0D+00 / real ( nt, kind = 8 ) cw = 0.0D+00 else if ( rule == 3 ) then ra(1) = 1.0D+00 rw(1) = 1.0D+00 call tvec_even ( nt, ta ) tw(1:4) = 0.125D+00 cw = 0.5D+00 else if ( rule == 4 ) then ra(1) = sqrt ( 2.0D+00 / 3.0D+00 ) rw(1) = 1.0D+00 call tvec_even ( nt, ta ) tw(1:nt) = 0.125D+00 cw = 0.25D+00 else if ( rule == 5 ) then a = 1.0D+00 b = sqrt ( 2.0D+00 ) / 2.0D+00 u = 1.0D+00 / 6.0D+00 v = 4.0D+00 / 6.0D+00 ra(1:nr) = (/ a, b /) rw(1:nr) = (/ u, v /) call tvec_even ( nt, ta ) tw(1:nt) = 1.0D+00 / real ( nt, kind = 8 ) cw = 4.0D+00 / 24.0D+00 else if ( rule == 6 ) then a = sqrt ( 3.0D+00 ) / 2.0D+00 b = sqrt ( ( 27.0D+00 - 3.0D+00 * sqrt ( 29.0D+00 ) ) / 52.0D+00 ) c = sqrt ( ( 27.0D+00 + 3.0D+00 * sqrt ( 29.0D+00 ) ) / 52.0D+00 ) u = 8.0D+00 / 27.0D+00 v = ( 551.0D+00 + 41.0D+00 * sqrt ( 29.0D+00 ) ) / 1566.0D+00 w = ( 551.0D+00 - 41.0D+00 * sqrt ( 29.0D+00 ) ) / 1566.0D+00 ra(1:nr) = (/ a, b, c /) rw(1:nr) = (/ u, v, w /) call tvec_even ( nt, ta ) tw(1:nt) = 1.0D+00 / real ( nt, kind = 8 ) cw = 0.0D+00 else if ( rule == 7 ) then a = sqrt ( ( 6.0D+00 - sqrt ( 6.0D+00 ) ) / 10.0D+00 ) b = sqrt ( ( 6.0D+00 + sqrt ( 6.0D+00 ) ) / 10.0D+00 ) u = ( 16.0D+00 + sqrt ( 6.0D+00 ) ) / 36.0D+00 v = ( 16.0D+00 - sqrt ( 6.0D+00 ) ) / 36.0D+00 ra(1:nr) = (/ a, b /) rw(1:nr) = (/ u, v /) call tvec_even ( nt, ta ) tw(1:nt) = 1.0D+00 / real ( nt, kind = 8 ) cw = 1.0D+00 / 9.0D+00 else if ( rule == 8 ) then call legendre_set ( nr, ra, rw ) a = -1.0D+00 b = +1.0D+00 c = 0.0D+00 d = +1.0D+00 call rule_adjust ( a, b, c, d, nr, ra, rw ) ra(1:nr) = sqrt ( ra(1:nr) ) call tvec_even ( nt, ta ) tw(1:nt) = 1.0D+00 / real ( nt, kind = 8 ) cw = 0.0D+00 else if ( rule == 9 ) then call legendre_set ( nr, ra, rw ) a = -1.0D+00 b = +1.0D+00 c = 0.0D+00 d = +1.0D+00 call rule_adjust ( a, b, c, d, nr, ra, rw ) ra(1:nr) = sqrt ( ra(1:nr) ) call tvec_even ( nt, ta ) tw(1:nt) = 1.0D+00 / real ( nt, kind = 8 ) cw = 0.0D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CIRCLE_RT_SET - Fatal error!' write ( *, '(a,i8)' ) ' There is no rule of index ', rule stop end if return end subroutine circle_rt_size ( rule, nr, nt, nc ) !*****************************************************************************80 ! !! CIRCLE_RT_SIZE sizes an R, THETA product quadrature rule in the unit circle. ! ! Discussion: ! ! For a given value of RULE, here are the number of points used at the ! center (NC), the number of points along the radial direction (NR) and ! the number of points along the theta direction (NT). The total number ! of points in the rule will be ! ! Total = NC + NR * NT. ! ! The user, when choosing RULE, must allocate enough space in the arrays ! RA, RW, TA and TW for the resulting values of NR and NT. ! ! RULE NC NR NT Total ! ---- -- -- -- ----- ! 1 1 0 0 1 ! 2 0 1 4 4 ! 3 1 1 4 5 ! 4 1 1 6 7 ! 5 1 2 4 9 ! 6 0 3 4 12 ! 7 1 2 10 21 ! 8 0 4 16 64 ! 9 0 5 20 120 ! ! The integral of F(X,Y) over the unit circle is approximated by ! ! Integral ( X*X + Y*Y <= 1 ) F(X,Y) dx dy ! = Integral ( 0 <= R <= 1, 0 <= T <= 2PI ) F(R*cos(T),R*sin(T)) r dr dt ! = approximately ! ZW * F(0,0) ! + sum ( 1 <= I <= NR ) Sum ( 1 <= J <= NT ) ! RW(I) * TW(J) * F ( R(I) * cos ( TA(J) ), R(I) * sin ( TA(J) ) ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 March 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964, ! ISBN: 0-486-61272-4, ! LC: QA47.A34. ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, integer ( kind = 4 ) RULE, the rule desired. ! ! Output, integer ( kind = 4 ) NR, the number of R abscissas. ! ! Output, integer ( kind = 4 ) NT, the number of Theta abscissas. ! ! Output, integer ( kind = 4 ) NC, the number of center abscissas (0 or 1). ! implicit none integer ( kind = 4 ) nc integer ( kind = 4 ) nr integer ( kind = 4 ) nt integer ( kind = 4 ) rule if ( rule == 1 ) then nr = 0 nt = 0 nc = 1 else if ( rule == 2 ) then nr = 1 nt = 4 nc = 0 else if ( rule == 3 ) then nr = 1 nt = 4 nc = 1 else if ( rule == 4 ) then nr = 1 nt = 6 nc = 1 else if ( rule == 5 ) then nr = 2 nt = 4 nc = 1 else if ( rule == 6 ) then nr = 3 nt = 4 nc = 0 else if ( rule == 7 ) then nr = 2 nt = 10 nc = 1 else if ( rule == 8 ) then nr = 4 nt = 16 nc = 0 else if ( rule == 9 ) then nr = 5 nt = 20 nc = 0 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CIRCLE_RT_SIZE - Fatal error!' write ( *, '(a,i8)' ) ' There is no rule of index ', rule stop end if return end subroutine circle_rt_sum ( func, center, radius, nr, ra, rw, nt, ta, tw, zw, & result ) !*****************************************************************************80 ! !! CIRCLE_RT_SUM applies an R, THETA product quadrature rule inside a circle. ! ! Integration region: ! ! (X-CENTER(1))^2 + (Y-CENTER(2))^2 <= RADIUS^2. ! ! Discussion: ! ! The product rule is assumed to be have the form: ! ! Integral_Approx = ZW * F(CENTER(1),CENTER(2)) + ! sum ( 1 <= IR <= NR ) Sum ( 1 <= IT <= NT ) ! RW(IR) * TW(IT) * F ( CENTER(1) + R(IR) * RADIUS * Cos ( TA(IT) ), ! CENTER(2) + R(IR) * RADIUS * Sin ( TA(IT) ) ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 15 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function of two variables which is to be integrated, ! of the form: ! function func ( x, y ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! ! Input, real ( kind = 8 ) CENTER(2), the center of the circle. ! ! Input, real ( kind = 8 ) RADIUS, the radius of the circle. ! ! Input, integer ( kind = 4 ) NR, the number of R abscissas. ! ! Input, real ( kind = 8 ) RA(NR), RW(NR), the R abscissas and weights. ! ! Input, integer ( kind = 4 ) NT, the number of Theta abscissas. ! ! Input, real ( kind = 8 ) TA(NT), TW(NT), the THETA abscissas and weights. ! ! Input, real ( kind = 8 ) ZW, the weight to use for the center. ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 integer ( kind = 4 ) nr integer ( kind = 4 ) nt real ( kind = 8 ) center(dim_num) real ( kind = 8 ) circle_area_2d real ( kind = 8 ), external :: func integer ( kind = 4 ) ir integer ( kind = 4 ) it real ( kind = 8 ) quad real ( kind = 8 ) ra(nr) real ( kind = 8 ) radius real ( kind = 8 ) rct real ( kind = 8 ) result real ( kind = 8 ) rst real ( kind = 8 ) rw(nr) real ( kind = 8 ) ta(nt) real ( kind = 8 ) tw(nt) real ( kind = 8 ) volume real ( kind = 8 ) x real ( kind = 8 ) y real ( kind = 8 ) zw quad = 0.0D+00 if ( zw /= 0.0D+00 ) then x = center(1) y = center(2) quad = quad + zw * func ( x, y ) end if do it = 1, nt rct = radius * cos ( ta(it) ) rst = radius * sin ( ta(it) ) do ir = 1, nr x = center(1) + ra(ir) * rct y = center(2) + ra(ir) * rst quad = quad + tw(it) * rw(ir) * func ( x, y ) end do end do volume = circle_area_2d ( radius ) result = quad * volume return end subroutine circle_sector ( func, center, radius, theta1, theta2, nr, result ) !*****************************************************************************80 ! !! CIRCLE_SECTOR approximates an integral in a circular sector. ! ! Discussion: ! ! A sector is contained within a circular arc and the lines joining each ! endpoint of the arc to the center of the circle. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied function of two ! variables which is to be integrated, of the form: ! function func ( x, y ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! ! Input, real ( kind = 8 ) CENTER(2), the center of the circle. ! ! Input, real ( kind = 8 ) RADIUS, the radius of the circle. ! ! Input, real ( kind = 8 ) THETA1, THETA2, the angles defining the sector. ! The sector is measured from THETA1 to THETA2. ! ! Input, integer ( kind = 4 ) NR, the number of radial values used in the approximation ! of the integral. NR must be at least 1. Higher values improve the ! accuracy of the integration, at the cost of more function evaluations. ! ! Output, real ( kind = 8 ) RESULT, the approximation to the integral. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 integer ( kind = 4 ) nr real ( kind = 8 ) a real ( kind = 8 ) area real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) center(dim_num) real ( kind = 8 ) circle_sector_area_2d real ( kind = 8 ) d real ( kind = 8 ), external :: func integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) nt real ( kind = 8 ) quad real ( kind = 8 ) ra(nr) real ( kind = 8 ) radius real ( kind = 8 ) result real ( kind = 8 ) rw(nr) real ( kind = 8 ) ta(4*nr) real ( kind = 8 ) theta1 real ( kind = 8 ) theta2 real ( kind = 8 ) tw(4*nr) real ( kind = 8 ) x real ( kind = 8 ) y ! ! Set the radial abscissas and weights. ! call legendre_set ( nr, ra, rw ) a = -1.0D+00 b = +1.0D+00 c = 0.0D+00 d = radius * radius call rule_adjust ( a, b, c, d, nr, ra, rw ) ra(1:nr) = sqrt ( ra(1:nr) ) rw(1:nr) = rw(1:nr) / radius / radius ! ! Pick angles evenly spaced between THETA1 and THETA2, but do not ! include the endpoints, and use a half interval for the first and last. ! nt = 4 * nr call tvec_even_bracket3 ( nt, theta1, theta2, ta ) tw(1:nt) = 1.0D+00 / real ( nt, kind = 8 ) ! ! Approximate the integral. ! quad = 0.0D+00 do i = 1, nr do j = 1, nt x = center(1) + ra(i) * cos ( ta(j) ) y = center(2) + ra(i) * sin ( ta(j) ) quad = quad + rw(i) * tw(j) * func ( x, y ) end do end do area = circle_sector_area_2d ( radius, theta1, theta2 ) result = quad * area return end function circle_sector_area_2d ( r, theta1, theta2 ) !*****************************************************************************80 ! !! CIRCLE_SECTOR_AREA_2D returns the area of a circular sector in 2D. ! ! Discussion: ! ! A sector is contained within a circular arc and the lines joining each ! endpoint of the arc to the center of the circle. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R, the radius of the circle. ! ! Input, real ( kind = 8 ) THETA1, THETA2, the angles of the rays ! that delimit the sector. ! ! Output, real ( kind = 8 ) CIRCLE_SECTOR_AREA_2D, the area of the sector. ! implicit none real ( kind = 8 ) circle_sector_area_2d real ( kind = 8 ) r real ( kind = 8 ) theta1 real ( kind = 8 ) theta2 circle_sector_area_2d = 0.50D+00 * r * r * ( theta2 - theta1 ) return end function circle_triangle_area_2d ( r, theta1, theta2 ) !*****************************************************************************80 ! !! CIRCLE_TRIANGLE_AREA_2D returns the area of a circle triangle in 2D. ! ! Discussion: ! ! A circle triangle is formed by drawing a circular arc, and considering ! the triangle formed by the endpoints of the arc plus the center of ! the circle. ! ! The normal situation is that 0 < ( THETA2 - THETA1 ) < PI. Outside ! this range, the triangle can actually have NEGATIVE area. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R, the radius of the circle. ! ! Input, real ( kind = 8 ) THETA1, THETA2, the angles of the rays that ! delimit the arc. ! ! Output, real ( kind = 8 ) CIRCLE_TRIANGLE_AREA_2D, the (signed) area ! of the triangle. ! implicit none real ( kind = 8 ) circle_triangle_area_2d real ( kind = 8 ) r real ( kind = 8 ) theta1 real ( kind = 8 ) theta2 circle_triangle_area_2d = 0.5D+00 * r * r * sin ( theta2 - theta1 ) return end subroutine circle_xy_set ( rule, order, xtab, ytab, weight ) !*****************************************************************************80 ! !! CIRCLE_XY_SET sets an XY quadrature rule inside the unit circle in 2D. ! ! Integration region: ! ! X*X + Y*Y <= 1.0. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 December 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964, ! ISBN: 0-486-61272-4, ! LC: QA47.A34. ! ! Frank Lether, ! A Generalized Product Rule for the Circle, ! SIAM Journal on Numerical Analysis, ! Volume 8, Number 2, June 1971, pages 249-253. ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, integer ( kind = 4 ) RULE, the rule desired. ! 1, 1 point 1-st degree; ! 2, 4 point 3-rd degree, Stroud S2:3-1; ! 3, 4 point 3-rd degree, Lether #1; ! 4, 4 point 3-rd degree, Stroud S2:3-2; ! 5, 5 point 3-rd degree; ! 6, 7 point 5-th degree; ! 7, 9 point 5-th degree; ! 8, 9 point 5-th degree, Lether #2; ! 9, 12 point 7-th degree; ! 10, 16 point 7-th degree, Lether #3; ! 11, 21 point 9-th degree, Stroud S2:9-3; ! 12, 25 point 9-th degree, Lether #4 (after correcting error); ! 13, 64 point 15-th degree Gauss product rule. ! ! Input, integer ( kind = 4 ) ORDER, the order of the desired rule. ! ! Output, real ( kind = 8 ) XTAB(ORDER), YTAB(ORDER), the abscissas of ! the rule. ! ! Output, real ( kind = 8 ) WEIGHT(ORDER), the ORDER weights of the rule. ! implicit none integer order real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) d real ( kind = 8 ) e real ( kind = 8 ) f real ( kind = 8 ) g real ( kind = 8 ) h integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) nr real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) r real ( kind = 8 ) ra(4) real ( kind = 8 ) rw(4) integer ( kind = 4 ) rule real ( kind = 8 ) s real ( kind = 8 ) u real ( kind = 8 ) v real ( kind = 8 ) w real ( kind = 8 ) w1 real ( kind = 8 ) w2 real ( kind = 8 ) w3 real ( kind = 8 ) w4 real ( kind = 8 ) w5 real ( kind = 8 ) w6 real ( kind = 8 ) w7 real ( kind = 8 ) w8 real ( kind = 8 ) w9 real ( kind = 8 ) weight(order) real ( kind = 8 ) xtab(order) real ( kind = 8 ) ytab(order) real ( kind = 8 ) z if ( rule == 1 ) then xtab(1) = 0.0D+00 ytab(1) = 0.0D+00 weight(1) = 1.0D+00 else if ( rule == 2 ) then a = 0.5D+00 b = 0.25D+00 z = 0.0D+00 xtab(1:4) = (/ a, -a, z, z /) ytab(1:4) = (/ z, z, a, -a /) weight(1:4) = (/ b, b, b, b /) else if ( rule == 3 ) then a = 0.5D+00 b = 0.25D+00 xtab(1:4) = (/ a, -a, -a, a /) ytab(1:4) = (/ a, a, -a, -a /) weight(1:4) = (/ b, b, b, b /) else if ( rule == 4 ) then a = sqrt ( 2.0D+00 ) / 2.0D+00 b = 0.25D+00 xtab(1:4) = (/ a, -a, -a, a /) ytab(1:4) = (/ a, a, -a, -a /) weight(1:4) = (/ b, b, b, b /) else if ( rule == 5 ) then a = 1.0D+00 b = 0.5D+00 c = 0.125D+00 z = 0.0D+00 xtab(1:5) = (/ z, a, z, -a, z /) ytab(1:5) = (/ z, z, a, z, -a /) weight(1:5) = (/ b, c, c, c, c /) else if ( rule == 6 ) then a = sqrt ( 2.0D+00 / 3.0D+00 ) b = sqrt ( 1.0D+00 / 6.0D+00 ) c = sqrt ( 2.0D+00 ) / 2.0D+00 d = 0.125D+00 e = 0.25D+00 z = 0.0D+00 xtab(1:7) = (/ z, a, -a, b, -b, b, -b /) ytab(1:7) = (/ z, z, z, c, c, -c, -c /) weight(1:7) = (/ e, d, d, d, d, d, d /) else if ( rule == 7 ) then a = 0.5D+00 b = 1.0D+00 c = 4.0D+00 / 24.0D+00 d = 1.0D+00 / 24.0D+00 z = 0.0D+00 xtab(1:9) = (/ z, b, -b, z, z, a, -a, -a, a /) ytab(1:9) = (/ z, z, z, b, -b, a, a, -a, -a /) weight(1:9) = (/ c, d, d, d, d, c, c, c, c /) else if ( rule == 8 ) then a = sqrt ( 2.0D+00 ) / 2.0D+00 b = sqrt ( 3.0D+00 / 5.0D+00 ) c = sqrt ( 3.0D+00 / 10.0D+00 ) w1 = 16.0D+00 / 72.0D+00 w2 = 8.0D+00 / 72.0D+00 w3 = 10.0D+00 / 72.0D+00 w4 = 5.0D+00 / 72.0D+00 z = 0.0D+00 xtab(1:9) = (/ z, a, -a, z, z, a, a, -a, -a /) ytab(1:9) = (/ z, z, z, b, -b, c, -c, c, -c /) weight(1:9) = (/ w1, w2, w2, w3, w3, w4, w4, w4, w4 /) else if ( rule == 9 ) then a = sqrt ( 3.0D+00 ) / 2.0D+00 b = sqrt ( ( 27.0D+00 - 3.0D+00 * sqrt ( 29.0D+00 ) ) / 104.0D+00 ) c = sqrt ( ( 27.0D+00 + 3.0D+00 * sqrt ( 29.0D+00 ) ) / 104.0D+00 ) u = 2.0D+00 / 27.0D+00 v = ( 551.0D+00 + 41.0D+00 * sqrt ( 29.0D+00 ) ) / 6264.0D+00 w = ( 551.0D+00 - 41.0D+00 * sqrt ( 29.0D+00 ) ) / 6264.0D+00 z = 0.0D+00 xtab(1:12) = (/ a, -a, z, z, b, -b, b, -b, c, c, -c, -c /) ytab(1:12) = (/ z, z, a, -a, b, b, -b, -b, c, -c, c, -c /) weight(1:12) = (/ u, u, u, u, v, v, v, v, w, w, w, w /) else if ( rule == 10 ) then a = sqrt ( ( 3.0D+00 - sqrt ( 5.0D+00 ) ) / 8.0D+00 ) b = sqrt ( ( 15.0D+00 + 3.0D+00 * sqrt ( 5.0D+00 ) & - 2.0D+00 * sqrt ( 30.0D+00 ) - 2.0D+00 * sqrt ( 6.0D+00 ) ) / 56.0D+00 ) c = sqrt ( ( 15.0D+00 + 3.0D+00 * sqrt ( 5.0D+00 ) & + 2.0D+00 * sqrt ( 30.0D+00 ) + 2.0D+00 * sqrt ( 6.0D+00 ) ) / 56.0D+00 ) d = sqrt ( ( 3.0D+00 + sqrt ( 5.0D+00 ) ) / 8.0D+00 ) e = sqrt ( ( 15.0D+00 - 3.0D+00 * sqrt ( 5.0D+00 ) & - 2.0D+00 * sqrt ( 30.0D+00 ) + 2.0D+00 * sqrt ( 6.0D+00 ) ) / 56.0D+00 ) f = sqrt ( ( 15.0D+00 - 3.0D+00 * sqrt ( 5.0D+00 ) & + 2.0D+00 * sqrt ( 30.0D+00 ) - 2.0D+00 * sqrt ( 6.0D+00 ) ) / 56.0D+00 ) w1 = ( 90.0D+00 + 5.0D+00 * sqrt ( 30.0D+00 ) + 18.0D+00 * sqrt ( 5.0D+00 ) & + 5.0D+00 * sqrt ( 6.0D+00 ) ) / 1440.0D+00 w2 = ( 90.0D+00 - 5.0D+00 * sqrt ( 30.0D+00 ) + 18.0D+00 * sqrt ( 5.0D+00 ) & - 5.0D+00 * sqrt ( 6.0D+00 ) ) / 1440.0D+00 w3 = ( 90.0D+00 + 5.0D+00 * sqrt ( 30.0D+00 ) - 18.0D+00 * sqrt ( 5.0D+00 ) & - 5.0D+00 * sqrt ( 6.0D+00 ) ) / 1440.0D+00 w4 = ( 90.0D+00 - 5.0D+00 * sqrt ( 30.0D+00 ) - 18.0D+00 * sqrt ( 5.0D+00 ) & + 5.0D+00 * sqrt ( 6.0D+00 ) ) / 1440.0D+00 xtab(1:order) = (/ a, a, -a, -a, a, a, -a, -a, d, d, -d, -d, & d, d, -d, -d /) ytab(1:order) = (/ b, -b, b, -b, c, -c, c, -c, e, -e, e, -e, & f, -f, f, -f /) weight(1:order) = (/ w1, w1, w1, w1, w2, w2, w2, w2, w3, w3, w3, w3, & w4, w4, w4, w4 /) else if ( rule == 11 ) then xtab(1) = 0.0D+00 ytab(1) = 0.0D+00 weight(1) = 1.0D+00 / 9.0D+00 weight(2:11) = ( 16.0D+00 + sqrt ( 6.0D+00 ) ) / 360.0D+00 weight(12:21) = ( 16.0D+00 - sqrt ( 6.0D+00 ) ) / 360.0D+00 r = sqrt ( ( 6.0D+00 - sqrt ( 6.0D+00 ) ) / 10.0D+00 ) do i = 1, 10 a = 2.0D+00 * pi * real ( i, kind = 8 ) / 10.0D+00 xtab(1+i) = r * cos ( a ) ytab(1+i) = r * sin ( a ) end do r = sqrt ( ( 6.0D+00 + sqrt ( 6.0D+00 ) ) / 10.0D+00 ) do i = 1, 10 a = 2.0D+00 * pi * real ( i, kind = 8 ) / 10.0D+00 xtab(11+i) = r * cos ( a ) ytab(11+i) = r * sin ( a ) end do ! ! There was apparently a misprint in the Lether paper. The quantity ! which here reads "322" was printed there as "332". ! else if ( rule == 12 ) then a = 0.5D+00 b = sqrt ( 3.0D+00 ) / 2.0D+00 c = sqrt ( ( 35.0D+00 + 2.0D+00 * sqrt ( 70.0D+00 ) ) / 252.0D+00 ) d = sqrt ( ( 35.0D+00 - 2.0D+00 * sqrt ( 70.0D+00 ) ) / 252.0D+00 ) e = sqrt ( ( 35.0D+00 + 2.0D+00 * sqrt ( 70.0D+00 ) ) / 84.0D+00 ) f = sqrt ( ( 35.0D+00 - 2.0D+00 * sqrt ( 70.0D+00 ) ) / 84.0D+00 ) g = sqrt ( ( 35.0D+00 + 2.0D+00 * sqrt ( 70.0D+00 ) ) / 63.0D+00 ) h = sqrt ( ( 35.0D+00 - 2.0D+00 * sqrt ( 70.0D+00 ) ) / 63.0D+00 ) w1 = 64.0D+00 / 675.0D+00 w2 = 16.0D+00 / 225.0D+00 w3 = 16.0D+00 / 675.0D+00 w4 = ( 322.0D+00 - 13.0D+00 * sqrt ( 70.0D+00 ) ) / 21600.0D+00 w5 = ( 322.0D+00 + 13.0D+00 * sqrt ( 70.0D+00 ) ) / 21600.0D+00 w6 = ( 322.0D+00 - 13.0D+00 * sqrt ( 70.0D+00 ) ) / 7200.0D+00 w7 = ( 322.0D+00 + 13.0D+00 * sqrt ( 70.0D+00 ) ) / 7200.0D+00 w8 = ( 322.0D+00 - 13.0D+00 * sqrt ( 70.0D+00 ) ) / 5400.0D+00 w9 = ( 322.0D+00 + 13.0D+00 * sqrt ( 70.0D+00 ) ) / 5400.0D+00 z = 0.0D+00 xtab(1:order) = (/ z, a, -a, b, -b, b, b, -b, -b, b, b, -b, -b, & a, a, -a, -a, a, a, -a, -a, z, z, z, z /) ytab(1:order) = (/ z, z, z, z, z, c, -c, c, -c, d, -d, d, -d, & e, -e, e, -e, f, -f, f, -f, g, -g, h, -h /) weight(1:order) = (/ w1, w2, w2, w3, w3, w4, w4, w4, w4, w5, w5, w5, w5, & w6, w6, w6, w6, w7, w7, w7, w7, w8, w8, w9, w9 /) else if ( rule == 13 ) then nr = 4 call legendre_set ( nr, ra, rw ) a = -1.0D+00 b = +1.0D+00 c = 0.0D+00 d = +1.0D+00 call rule_adjust ( a, b, c, d, nr, ra, rw ) ra(1:nr) = sqrt ( ra(1:nr) ) i = 0 do j = 1, 16 c = cos ( pi * real ( j, kind = 8 ) / 8.0D+00 ) s = sin ( pi * real ( j, kind = 8 ) / 8.0D+00 ) do k = 1, 4 i = i + 1 xtab(i) = c * ra(k) ytab(i) = s * ra(k) weight(i) = rw(k) / 16.0D+00 end do end do else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CIRCLE_XY_SET - Fatal error!' write ( *, '(a,i8)' ) ' There is no rule of index ', rule stop end if return end subroutine circle_xy_size ( rule, order ) !*****************************************************************************80 ! !! CIRCLE_XY_SIZE sizes an XY quadrature rule inside the unit circle in 2D. ! ! Integration region: ! ! X*X + Y*Y <= 1.0. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 March 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964, ! ISBN: 0-486-61272-4, ! LC: QA47.A34. ! ! Frank Lether, ! A Generalized Product Rule for the Circle, ! SIAM Journal on Numerical Analysis, ! Volume 8, Number 2, June 1971, pages 249-253. ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, integer ( kind = 4 ) RULE, the rule desired. ! 1, 1 point 1-st degree; ! 2, 4 point 3-rd degree, Stroud S2:3-1; ! 3, 4 point 3-rd degree, Lether #1; ! 4, 4 point 3-rd degree, Stroud S2:3-2; ! 5, 5 point 3-rd degree; ! 6, 7 point 5-th degree; ! 7, 9 point 5-th degree; ! 8, 9 point 5-th degree, Lether #2; ! 9, 12 point 7-th degree; ! 10, 16 point 7-th degree, Lether #3; ! 11, 21 point 9-th degree, Stroud S2:9-3; ! 12, 25 point 9-th degree, Lether #4 (after correcting error); ! 13, 64 point 15-th degree Gauss product rule. ! ! Output, integer ( kind = 4 ) ORDER, the order of the desired rule. ! implicit none integer ( kind = 4 ) order integer ( kind = 4 ) rule if ( rule == 1 ) then order = 1 else if ( rule == 2 ) then order = 4 else if ( rule == 3 ) then order = 4 else if ( rule == 4 ) then order = 4 else if ( rule == 5 ) then order = 5 else if ( rule == 6 ) then order = 7 else if ( rule == 7 ) then order = 9 else if ( rule == 8 ) then order = 9 else if ( rule == 9 ) then order = 12 else if ( rule == 10 ) then order = 16 else if ( rule == 11 ) then order = 21 else if ( rule == 12 ) then order = 25 else if ( rule == 13 ) then order = 64 else order = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CIRCLE_XY_SIZE - Fatal error!' write ( *, '(a,i8)' ) ' There is no rule of index ', rule stop end if return end subroutine circle_xy_sum ( func, center, r, order, xtab, ytab, weight, & result ) !*****************************************************************************80 ! !! CIRCLE_XY_SUM applies an XY quadrature rule inside a circle in 2D. ! ! Integration region: ! ! (X-CENTER(1))^2 + (Y-CENTER(2))^2 <= R * R. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function of two variables which is to be integrated, ! of the form: ! function func ( x, y ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! ! Input, real ( kind = 8 ) CENTER(2), the coordinates of the center of ! the circle. ! ! Input, real ( kind = 8 ) R, the radius of the circle. ! ! Input, integer ( kind = 4 ) ORDER, the order of the rule. The rule is ! assumed to be defined on the unit circle. ! ! Input, real ( kind = 8 ) XTAB(ORDER), YTAB(ORDER), the XY ! coordinates of the abscissas of the quadrature rule for the unit circle. ! ! Input, real ( kind = 8 ) WEIGHT(ORDER), the weights of the rule. ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 integer ( kind = 4 ) order real ( kind = 8 ) center(dim_num) real ( kind = 8 ) circle_area_2d real ( kind = 8 ), external :: func integer ( kind = 4 ) i real ( kind = 8 ) quad real ( kind = 8 ) r real ( kind = 8 ) result real ( kind = 8 ) volume real ( kind = 8 ) weight(order) real ( kind = 8 ) x real ( kind = 8 ) xtab(order) real ( kind = 8 ) y real ( kind = 8 ) ytab(order) quad = 0.0D+00 do i = 1, order x = center(1) + r * xtab(i) y = center(2) + r * ytab(i) quad = quad + weight(i) * func ( x, y ) end do volume = circle_area_2d ( r ) result = quad * volume return end subroutine cone_unit_3d ( func, result ) !*****************************************************************************80 ! !! CONE_UNIT_3D approximates an integral inside the unit cone in 3D. ! ! Integration Region: ! ! X*X + Y*Y <= 1 - Z ! ! and ! ! 0 <= Z <= 1. ! ! Discussion: ! ! An 48 point degree 7 formula, Stroud CN:S2:7-1, is used. ! ! (There is a typographical error in the S2:7-1 formula for B3.) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 18 November 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied function which ! evaluates F(X,Y,Z), of the form ! function func ( x, y, z ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! real ( kind = 8 ) z ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) cone_volume_3d real ( kind = 8 ), external :: func real ( kind = 8 ) h integer ( kind = 4 ) i real ( kind = 8 ) quad real ( kind = 8 ) r real ( kind = 8 ) result real ( kind = 8 ), save, dimension ( 4 ) :: u = & (/ 0.04850054945D+00, 0.2386007376D+00, & 0.5170472951D+00, 0.7958514179D+00 /) real ( kind = 8 ) volume real ( kind = 8 ), save, dimension ( 4 ) :: w1 = & (/ 0.1108884156D+00, 0.1434587878D+00, & 0.06863388717D+00, 0.01035224075D+00 /) real ( kind = 8 ) w2(3) real ( kind = 8 ) x real ( kind = 8 ) y real ( kind = 8 ) z a = sqrt ( 3.0D+00 ) / 2.0D+00 b = sqrt ( ( 27.0D+00 - 3.0D+00 * sqrt ( 29.0D+00 ) ) / 104.0D+00 ) c = sqrt ( ( 27.0D+00 + 3.0D+00 * sqrt ( 29.0D+00 ) ) / 104.0D+00 ) w2(1:3) = 3.0D+00 * (/ & 2.0D+00 / 27.0D+00, & ( 551.0D+00 + 4.0D+00 * sqrt ( 29.0D+00 ) ) / 6264.0D+00, & ( 551.0D+00 - 4.0D+00 * sqrt ( 29.0D+00 ) ) / 6264.0D+00 /) quad = 0.0D+00 do i = 1, 4 x = a * ( 1.0D+00 - u(i) ) y = 0.0D+00 z = u(i) quad = quad + w1(i) * w2(1) * func ( x, y, z ) x = -a * ( 1.0D+00 - u(i) ) y = 0.0D+00 z = u(i) quad = quad + w1(i) * w2(1) * func ( x, y, z ) x = 0.0D+00 y = a * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(1) * func ( x, y, z ) x = 0.0D+00 y = -a * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(1) * func ( x, y, z ) end do do i = 1, 4 x = b * ( 1.0D+00 - u(i) ) y = b * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(2) * func ( x, y, z ) x = -b * ( 1.0D+00 - u(i) ) y = b * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(2) * func ( x, y, z ) x = -b * ( 1.0D+00 - u(i) ) y = -b * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(2) * func ( x, y, z ) x = b * ( 1.0D+00 - u(i) ) y = -b * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(2) * func ( x, y, z ) x = c * ( 1.0D+00 - u(i) ) y = c * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(3) * func ( x, y, z ) x = -c * ( 1.0D+00 - u(i) ) y = c * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(3) * func ( x, y, z ) x = -c * ( 1.0D+00 - u(i) ) y = -c * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(3) * func ( x, y, z ) x = c * ( 1.0D+00 - u(i) ) y = -c * ( 1.0D+00 - u(i) ) z = u(i) quad = quad + w1(i) * w2(3) * func ( x, y, z ) end do r = 1.0D+00 h = 1.0D+00 volume = cone_volume_3d ( r, h ) result = quad * volume return end function cone_volume_3d ( r, h ) !*****************************************************************************80 ! !! CONE_VOLUME_3D returns the volume of a cone in 3D. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R, the radius of the base of the cone. ! ! Input, real ( kind = 8 ) H, the height of the cone. ! ! Output, real ( kind = 8 ) CONE_VOLUME_3D, the volume of the cone. ! implicit none real ( kind = 8 ) cone_volume_3d real ( kind = 8 ) h real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) r cone_volume_3d = ( pi / 3.0D+00 ) * h * r * r return end subroutine cube_shell_nd ( func, n, r1, r2, result ) ! !*****************************************************************************80 ! !! CUBE_SHELL_ND approximates an integral inside a cubic shell in N dimensions. ! ! Integration region: ! ! R1 <= abs ( X(1:N) ) <= R2 ! ! Discussion: ! ! An N*2**N point third degree formula is used, Stroud number CNSHELL:3-4. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 31 May 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates F at the N-vector X, of the form ! function func ( n, x ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) func ! real ( kind = 8 ) x(n) ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Input, real ( kind = 8 ) R1, R2, the inner and outer radii of the cubical ! shell. The outer cube is of side 2*R2, the inner, missing cube of side ! 2*R1. ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) cube_shell_volume_nd logical done real ( kind = 8 ), external :: func integer ( kind = 4 ) i real ( kind = 8 ) quad real ( kind = 8 ) r1 real ( kind = 8 ) r2 real ( kind = 8 ) rmax real ( kind = 8 ) rmin real ( kind = 8 ) result real ( kind = 8 ) u real ( kind = 8 ) v real ( kind = 8 ) volume real ( kind = 8 ) x(n) if ( r1 == r2 ) then result = 0.0D+00 return end if rmax = max ( r1, r2 ) rmin = min ( r1, r2 ) u = sqrt ( real ( n, kind = 8 ) * ( rmax**(n+2) - rmin**(n+2) ) & / ( real ( n + 2, kind = 8 ) * ( rmax**n - rmin**n ) ) ) v = u / sqrt ( 3.0D+00 ) quad = 0.0D+00 do i = 1, n x(1:n) = v x(i) = u do quad = quad + func ( n, x ) call r8vec_mirror_next ( n, x, done ) if ( done ) then exit end if end do end do quad = quad / real ( n * 2**n, kind = 8 ) volume = cube_shell_volume_nd ( n, r1, r2 ) result = quad * volume return end function cube_shell_volume_nd ( n, r1, r2 ) !*****************************************************************************80 ! !! CUBE_SHELL_VOLUME_ND computes the volume of a cubic shell in ND. ! ! Integration region: ! ! R1 <= abs ( X(1:N) ) <= R2 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 20 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Input, real ( kind = 8 ) R1, R2, the inner and outer radii of the cubic ! shell. The outer cube is of side 2*R2, the inner, missing cube of side ! 2*R1. ! ! Output, real ( kind = 8 ) CUBE_SHELL_VOLUME_ND, the volume of the cubic ! shell. ! implicit none real ( kind = 8 ) cube_shell_volume_nd integer ( kind = 4 ) n real ( kind = 8 ) r1 real ( kind = 8 ) r2 cube_shell_volume_nd = ( r2**n - r1**n ) * 2**n return end subroutine cube_unit_3d ( func, result ) !*****************************************************************************80 ! !! CUBE_UNIT_3D approximates an integral inside the unit cube in 3D. ! ! Integration region: ! ! -1 <= X <= 1, ! and ! -1 <= Y <= 1, ! and ! -1 <= Z <= 1. ! ! Discussion: ! ! An 8 point third degree formula is used. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 15 August 1998 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates F(X,Y,Z), of the form ! function func ( x, y, z ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! real ( kind = 8 ) z ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none real ( kind = 8 ) cube_unit_volume_nd real ( kind = 8 ), external :: func real ( kind = 8 ) quad real ( kind = 8 ) result real ( kind = 8 ) s real ( kind = 8 ) volume real ( kind = 8 ) w real ( kind = 8 ) x real ( kind = 8 ) y real ( kind = 8 ) z s = 1.0D+00 / sqrt ( 3.0D+00 ) w = 1.0D+00 / 8.0D+00 x = s y = s z = s quad = w * ( & func ( x, y, z ) + func ( x, y, -z ) & + func ( x, -y, z ) + func ( x, -y, -z ) & + func ( -x, y, z ) + func ( -x, y, -z ) & + func ( -x, -y, z ) + func ( -x, -y, -z ) ) volume = cube_unit_volume_nd ( 3 ) result = quad * volume return end subroutine cube_unit_nd ( func, qa, qb, n, k ) !*****************************************************************************80 ! !! CUBE_UNIT_ND approximates an integral inside the unit cube in ND. ! ! Integration region: ! ! -1 <= X(1:N) <= 1 ! ! Discussion: ! ! A K**N point product formula is used. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 March 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! James Lyness, BJJ McHugh, ! Integration Over Multidimensional Hypercubes, ! A Progressive Procedure, ! The Computer Journal, ! Volume 6, 1963, pages 264-270. ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function which evaluates the function, of the form ! function func ( n, x ) ! integer ( kind = 4 ) n ! real ( kind = 8 ) func ! real ( kind = 8 ) x(n) ! ! Output, real ( kind = 8 ) QA(K), QB(K), two sets of estimates for ! the integral. The QB entries are obtained from the ! QA entries by Richardson extrapolation, and QB(K) is ! the best estimate for the integral. ! ! Input, integer ( kind = 4 ) N, the dimension of the cube. ! ! Input, integer ( kind = 4 ) K, the highest order of integration, and the order ! of Richardson extrapolation. K can be no greater than 10. ! implicit none integer ( kind = 4 ), parameter :: kmax = 10 integer ( kind = 4 ) k integer ( kind = 4 ) n real ( kind = 8 ), external :: func real ( kind = 8 ), save, dimension ( kmax, kmax ) :: g integer ( kind = 4 ) i real ( kind = 8 ) qa(k) real ( kind = 8 ) qb(k) g(1:kmax,1:kmax) = 0.0D+00 g( 1, 1) = 1.0D+00 g( 2, 1) = -0.3333333333333D+00 g( 2, 2) = 0.1333333333333D+01 g( 3, 1) = 0.4166666666667D-01 g( 3, 2) = -0.1066666666667D+01 g( 3, 3) = 0.2025000000000D+01 g( 4, 1) = -0.2777777777778D-02 g( 4, 2) = 0.3555555555556D+00 g( 4, 3) = -0.2603571428571D+01 g( 4, 4) = 0.3250793650794D+01 g( 5, 1) = 0.1157407407407D-03 g( 5, 2) = -0.6772486772487D-01 g( 5, 3) = 0.1464508928571D+01 g( 5, 4) = -0.5779188712522D+01 g( 5, 5) = 0.5382288910935D+01 g( 6, 1) = -0.3306878306878D-05 g( 6, 2) = 0.8465608465608D-02 g( 6, 3) = -0.4881696428571D+00 g( 6, 4) = 0.4623350970018D+01 g( 6, 5) = -0.1223247479758D+02 g( 6, 6) = 0.9088831168831D+01 g( 7, 1) = 0.6889329805996D-07 g( 7, 2) = -0.7524985302763D-03 g( 7, 3) = 0.1098381696429D+00 g( 7, 4) = -0.2241624712736D+01 g( 7, 5) = 0.1274216124748D+02 g( 7, 6) = -0.2516907092907D+02 g( 7, 7) = 0.1555944865432D+02 g( 8, 1) = -0.1093544413650D-08 g( 8, 2) = 0.5016656868509D-04 g( 8, 3) = -0.1797351866883D-01 g( 8, 4) = 0.7472082375786D+00 g( 8, 5) = -0.8168052081717D+01 g( 8, 6) = 0.3236023405166D+02 g( 8, 7) = -0.5082753227079D+02 g( 8, 8) = 0.2690606541646D+02 g( 9, 1) = 0.1366930517063D-10 g( 9, 2) = -0.2606055516108D-05 g( 9, 3) = 0.2246689833604D-02 g( 9, 4) = -0.1839281815578D+00 g( 9, 5) = 0.3646451822195D+01 g( 9, 6) = -0.2588818724133D+02 g( 9, 7) = 0.7782965878964D+02 g( 9, 8) = -0.1012934227443D+03 g( 9, 9) = 0.4688718347156D+02 g(10, 1) = -0.1380737896023D-12 g(10, 2) = 0.1085856465045D-06 g(10, 3) = -0.2222000934334D-03 g(10, 4) = 0.3503393934435D-01 g(10, 5) = -0.1215483940732D+01 g(10, 6) = 0.1456210532325D+02 g(10, 7) = -0.7477751530769D+02 g(10, 8) = 0.1800771959898D+03 g(10, 9) = -0.1998874663788D+03 g(10,10) = 0.8220635246624D+02 if ( kmax < k ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CUBE_UNIT_ND - Fatal error!' write ( *, '(a,i8)' ) ' K must be no greater than KMAX = ', kmax write ( *, '(a,i8)' ) ' but the input K is ', k stop end if do i = 1, k call qmdpt ( func, n, i, qa(i) ) end do qb(1) = qa(1) do i = 2, k qb(i) = dot_product ( g(i,1:i), qa(1:i) ) end do return end function cube_unit_volume_nd ( n ) !*****************************************************************************80 ! !! CUBE_UNIT_VOLUME_ND returns the volume of the unit cube in ND. ! ! Integration region: ! ! -1 <= X(1:N) <= 1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the dimension of the space. ! ! Output, real ( kind = 8 ) CUBE_UNIT_VOLUME_ND, the volume of the unit ! cube in ND. ! implicit none real ( kind = 8 ) cube_unit_volume_nd integer ( kind = 4 ) n cube_unit_volume_nd = 2.0D+00**n return end function ellipse_area_2d ( r1, r2 ) !*****************************************************************************80 ! !! ELLIPSE_AREA_2D returns the area of an ellipse in 2D. ! ! Integration region: ! ! ( ( X - CENTER(1) ) / R1 )^2 + ( ( Y - CENTER(2) ) / R2 )^2 <= 1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R1, R2, the major and minor semi-axes. ! ! Output, real ( kind = 8 ) ELLIPSE_AREA_2D, the area of the ellipse. ! implicit none real ( kind = 8 ) ellipse_area_2d real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) r1 real ( kind = 8 ) r2 ellipse_area_2d = pi * r1 * r2 return end function ellipse_circumference_2d ( r1, r2 ) !*****************************************************************************80 ! !! ELLIPSE_CIRCUMFERENCE_2D returns the circumference of an ellipse in 2D. ! ! Discussion: ! ! There is no closed formula for the circumference of an ellipse. ! ! Defining the eccentricity by ! ! E = sqrt ( 1 - ( r2 / r1 )^2 ) ! ! where R1 and R2 are the major and minor axes, then ! ! circumference ! = 4 * R1 * E(K,2*PI) ! = R1 * Integral ( 0 <= T <= 2*PI ) sqrt ( 1 - E * E * sin^2 ( T ) ) dT ! ! This integral can be approximated by the Gauss-Kummer formula. ! ! Integration region: ! ! ( ( X - CENTER(1) ) / R1 )^2 + ( ( Y - CENTER(2) ) / R2 )^2 <= 1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 21 May 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! John Harris, Horst Stocker, ! Handbook of Mathematics and Computational Science, ! Springer, 1998, ! ISBN: 0-387-94746-9, ! LC: QA40.S76. ! ! Parameters: ! ! Input, real ( kind = 8 ) R1, R2, the major and minor semi-axes. ! ! Output, real ( kind = 8 ) ELLIPSE_CIRCUMFERENCE_2D, the ! circumference of the ellipse. ! implicit none real ( kind = 8 ) ellipse_circumference_2d real ( kind = 8 ) e integer ( kind = 4 ) i real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) r1 real ( kind = 8 ) r2 real ( kind = 8 ) term real ( kind = 8 ) value if ( r1 == r2 ) then ellipse_circumference_2d = 2.0D+00 * pi * r1 return end if ! ! Compute the eccentricity of the ellipse. ! e = sqrt ( 1.0D+00 - ( min ( r1, r2 ) / max ( r1, r2 ) )**2 ) value = 1.0D+00 term = value i = 0 do i = i + 1 term = term * ( 2 * i - 3 ) * ( 2 * i - 1 ) * e * e & / real ( 2 * 2 * i * i, kind = 8 ) if ( abs ( term ) <= epsilon ( value ) * ( abs ( value ) + 1.0D+00 ) ) then exit end if value = value + term end do ellipse_circumference_2d = 2.0D+00 * pi * max ( r1, r2 ) * value return end function ellipse_eccentricity_2d ( r1, r2 ) !*****************************************************************************80 ! !! ELLIPSE_ECCENTRICITY_2D returns the eccentricity of an ellipse in 2D. ! ! Integration region: ! ! ( ( X - CENTER(1) ) / R1 )^2 + ( ( Y - CENTER(2) ) / R2 )^2 <= 1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R1, R2, the major and minor semi-axes. ! ! Output, real ( kind = 8 ) ELLIPSE_ECCENTRICITY_2D, the eccentricity ! of the ellipse. ! implicit none real ( kind = 8 ) ellipse_eccentricity_2d real ( kind = 8 ) major real ( kind = 8 ) minor real ( kind = 8 ) r1 real ( kind = 8 ) r2 minor = min ( abs ( r1 ), abs ( r2 ) ) major = max ( abs ( r1 ), abs ( r2 ) ) if ( major == 0.0D+00 ) then ellipse_eccentricity_2d = -huge ( r1 ) return end if ellipse_eccentricity_2d = sqrt ( 1.0D+00 - ( minor / major )**2 ) return end function ellipsoid_volume_3d ( r1, r2, r3 ) !*****************************************************************************80 ! !! ELLIPSOID_VOLUME_3D returns the volume of an ellipsoid in 3d. ! ! Discussion: ! ! This is not a general ellipsoid, but one for which each of the ! axes lies along a coordinate axis. ! ! Integration region: ! ! ( ( X - CENTER(1) ) / R1 )^2 ! + ( ( Y - CENTER(2) ) / R2 )^2 ! + ( ( Z - CENTER(3) ) / R3 )^2 <= 1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 May 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R1, R2, R3, the semi-axes of the ellipsoid. ! ! Output, real ( kind = 8 ) ELLIPSOID_VOLUME_3D, the volume of the ellipsoid. ! implicit none real ( kind = 8 ) ellipsoid_volume_3d real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) r1 real ( kind = 8 ) r2 real ( kind = 8 ) r3 ellipsoid_volume_3d = ( 4.0D+00 / 3.0D+00 ) * pi * r1 * r2 * r3 return end function hexagon_area_2d ( r ) !*****************************************************************************80 ! !! HEXAGON_AREA_2D returns the area of a regular hexagon in 2D. ! ! Discussion: ! ! The formula for the area only requires the radius, and does ! not depend on the location of the center, or the orientation ! of the hexagon. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) R, the radius of the hexagon. ! ! Output, real ( kind = 8 ) HEXAGON_AREA_2D, the area of the hexagon. ! implicit none real ( kind = 8 ) hexagon_area_2d real ( kind = 8 ) hexagon_unit_area_2d real ( kind = 8 ) r hexagon_area_2d = r * r * hexagon_unit_area_2d ( ) return end subroutine hexagon_sum ( func, center, r, order, xtab, ytab, weight, & result ) !*****************************************************************************80 ! !! HEXAGON_SUM applies a quadrature rule inside a hexagon in 2D. ! ! Integration region: ! ! The definition is given in terms of THETA, the angle in degrees of the ! vector (X-CENTER(1),Y-CENTER(2)). The following six conditions apply, ! respectively, between the bracketing values of THETA of 0, 60, 120, ! 180, 240, 300, and 360. ! ! 0 <= Y-CENTER(2) <= -SQRT(3) * (X-CENTER(1)) + R * SQRT(3) ! 0 <= Y-CENTER(2) <= R * SQRT(3)/2 ! 0 <= Y-CENTER(2) <= SQRT(3) * (X-CENTER(1)) + R * SQRT(3) ! -SQRT(3) * (X-CENTER(1)) - R * SQRT(3) <= Y-CENTER(2) <= 0 ! - R * SQRT(3)/2 <= Y-CENTER(2) <= 0 ! SQRT(3) * (X-CENTER(1)) - R * SQRT(3) <= Y-CENTER(2) <= 0 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, external FUNC, the name of the user supplied ! function of two variables which is to be integrated, ! of the form: ! function func ( x, y ) ! real ( kind = 8 ) func ! real ( kind = 8 ) x ! real ( kind = 8 ) y ! ! Input, real ( kind = 8 ) CENTER(2), the center of the hexagon. ! ! Input, real ( kind = 8 ) R, the radius of the hexagon. ! ! Input, integer ( kind = 4 ) ORDER, the order of the rule. ! ! Input, real ( kind = 8 ) XTAB(ORDER), YTAB(ORDER), the abscissas. ! ! Input, real ( kind = 8 ) WEIGHT(ORDER), the weights of the rule. ! ! Output, real ( kind = 8 ) RESULT, the approximate integral of the function. ! implicit none integer ( kind = 4 ), parameter :: dim_num = 2 integer ( kind = 4 ) order real ( kind = 8 ) center(dim_num) real ( kind = 8 ), external :: func real ( kind = 8 ) hexagon_area_2d integer ( kind = 4 ) i real ( kind = 8 ) quad real ( kind = 8 ) r real ( kind = 8 ) result real ( kind = 8 ) volume real ( kind = 8 ) weight(order) real ( kind = 8 ) x real ( kind = 8 ) xtab(order) real ( kind = 8 ) y real ( kind = 8 ) ytab(order) quad = 0.0D+00 do i = 1, order x = center(1) + r * xtab(i) y = center(2) + r * ytab(i) quad = quad + weight(i) * func ( x, y ) end do volume = hexagon_area_2d ( r ) result = quad * volume return end function hexagon_unit_area_2d ( ) !*****************************************************************************80 ! !! HEXAGON_UNIT_AREA_2D returns the area of the unit regular hexagon in 2D. ! ! Integration region: ! ! The definition is given in terms of THETA, the angle in degrees of the ! vector (X,Y). The following six conditions apply, respectively, ! between the bracketing values of THETA of 0, 60, 120, 180, 240, ! 300, and 360. ! ! 0 <= Y <= -SQRT(3) * X + SQRT(3) ! 0 <= Y <= SQRT(3)/2 ! 0 <= Y <= SQRT(3) * X + SQRT(3) ! - SQRT(3) * X - SQRT(3) <= Y <= 0 ! - SQRT(3)/2 <= Y <= 0 ! SQRT(3) * X - SQRT(3) <= Y <= 0 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) HEXAGON_UNIT_AREA_2D, the area of the hexagon. ! implicit none real ( kind = 8 ) hexagon_unit_area_2d hexagon_unit_area_2d = 3.0D+00 * sqrt ( 3.0D+00 ) / 2.0D+00 return end subroutine hexagon_unit_set ( rule, order, xtab, ytab, weight ) !*****************************************************************************80 ! !! HEXAGON_UNIT_SET sets a quadrature rule inside the unit hexagon in 2D. ! ! Integration region: ! ! The definition is given in terms of THETA, the angle in degrees of the ! vector (X,Y). The following six conditions apply, respectively, ! between the bracketing values of THETA of 0, 60, 120, 180, 240, ! 300, and 360. ! ! 0 <= Y <= -SQRT(3) * X + SQRT(3) ! 0 <= Y <= SQRT(3)/2 ! 0 <= Y <= SQRT(3) * X + SQRT(3) ! -SQRT(3) * X - SQRT(3) <= Y <= 0 ! - SQRT(3)/2 <= Y <= 0 ! SQRT(3) * X - SQRT(3) <= Y <= 0 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 09 December 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964, ! ISBN: 0-486-61272-4, ! LC: QA47.A34. ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, integer ( kind = 4 ) RULE, the rule desired. ! 1, 1 point, degree 1; ! 2, 4 points, degree 3; ! 3, 7 points, degree 3; ! 4, 7 points, degree 5; ! ! Input, integer ( kind = 4 ) ORDER, the order of the desired rule. ! ! Output, real ( kind = 8 ) XTAB(*), YTAB(*), the abscissas of the rule. ! ! Output, real ( kind = 8 ) WEIGHT(*), the ORDER weights of the rule. ! implicit none integer ( kind = 4 ) order real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) d real ( kind = 8 ) e integer ( kind = 4 ) rule real ( kind = 8 ) weight(order) real ( kind = 8 ) xtab(order) real ( kind = 8 ) ytab(order) real ( kind = 8 ) z if ( rule == 1 ) then xtab(1) = 0.0D+00 ytab(1) = 0.0D+00 weight(1) = 1.0D+00 ! ! Stroud rule H2:3-1. ! else if ( rule == 2 ) then a = sqrt ( 5.0D+00 / 12.0D+00 ) b = 1.0D+00 / 4.0D+00 z = 0.0D+00 xtab(1:order) = (/ a, -a, z, z /) ytab(1:order) = (/ z, z, a, -a /) weight(1:order) = (/ b, b, b, b /) ! ! Stroud rule H2:3-2. ! else if ( rule == 3 ) then a = sqrt ( 3.0D+00 ) / 2.0D+00 b = 0.5D+00 c = 1.0D+00 d = 5.0D+00 / 72.0D+00 e = 42.0D+00 / 72.0D+00 z = 0.0D+00 xtab(1:order) = (/ z, c, -c, b, -b, b, -b /) ytab(1:order) = (/ z, z, z, a, a, -a, -a /) weight(1:order) = (/ e, d, d, d, d, d, d /) ! ! Stroud rule H2:5-1. ! else if ( rule == 4 ) then a = sqrt ( 14.0D+00 ) / 5.0D+00 b = sqrt ( 14.0D+00 ) / 10.0D+00 c = sqrt ( 42.0D+00 ) / 10.0D+00 d = 125.0D+00 / 1008.0D+00 e = 258.0D+00 / 1008.0D+00 z = 0.0D+00 xtab(1:order) = (/ z, a, -a, b, -b, b, -b /) ytab(1:order) = (/ z, z, z, c, c, -c, -c /) weight(1:order) = (/ e, d, d, d, d, d, d /) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'HEXAGON_UNIT_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal input value of RULE = ', rule stop end if return end subroutine hexagon_unit_size ( rule, order ) !*****************************************************************************80 ! !! HEXAGON_UNIT_SIZE sizes a quadrature rule inside the unit hexagon in 2D. ! ! Integration region: ! ! The definition is given in terms of THETA, the angle in degrees of the ! vector (X,Y). The following six conditions apply, respectively, ! between the bracketing values of THETA of 0, 60, 120, 180, 240, ! 300, and 360. ! ! 0 <= Y <= -SQRT(3) * X + SQRT(3) ! 0 <= Y <= SQRT(3)/2 ! 0 <= Y <= SQRT(3) * X + SQRT(3) ! -SQRT(3) * X - SQRT(3) <= Y <= 0 ! - SQRT(3)/2 <= Y <= 0 ! SQRT(3) * X - SQRT(3) <= Y <= 0 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 13 March 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964, ! ISBN: 0-486-61272-4, ! LC: QA47.A34. ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, integer ( kind = 4 ) RULE, the rule desired. ! 1, 1 point, degree 1; ! 2, 4 points, degree 3; ! 3, 7 points, degree 3; ! 4, 7 points, degree 5; ! ! Output, integer ( kind = 4 ) ORDER, the order of the desired rule. ! If RULE is not legal, then ORDER is returned as -1. ! implicit none integer ( kind = 4 ) order integer ( kind = 4 ) rule if ( rule == 1 ) then order = 1 ! ! Stroud rule H2:3-1. ! else if ( rule == 2 ) then order = 4 ! ! Stroud rule H2:3-2. ! else if ( rule == 3 ) then order = 7 ! ! Stroud rule H2:5-1. ! else if ( rule == 4 ) then order = 7 else order = -1 end if return end function i4_factorial ( n ) !*****************************************************************************80 ! !! I4_FACTORIAL computes N! (for small values of N). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 08 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the argument. ! ! Output, integer ( kind = 4 ) I4_FACTORIAL, the value of N!. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) n integer ( kind = 4 ) i4_factorial i4_factorial = 1 do i = 1, n i4_factorial = i4_factorial * i end do return end function i4_factorial2 ( n ) !*****************************************************************************80 ! !! I4_FACTORIAL2 computes the factorial N!! ! ! Formula: ! ! FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 ) (N even) ! = Product ( N * (N-2) * (N-4) * ... * 1 ) (N odd) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 December 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the argument of the double factorial function. ! If N is less than 1, I4_FACTORIAL2 is returned as 1. ! ! Output, integer ( kind = 4 ) I4_FACTORIAL2, the value of N!!. ! implicit none integer ( kind = 4 ) i4_factorial2 integer ( kind = 4 ) n integer ( kind = 4 ) n_copy if ( n < 1 ) then i4_factorial2 = 1 return end if n_copy = n i4_factorial2 = 1 do while ( 1 < n_copy ) i4_factorial2 = i4_factorial2 * n_copy n_copy = n_copy - 2 end do return end subroutine ksub_next2 ( n, k, iarray, in, iout ) !*****************************************************************************80 ! !! KSUB_NEXT2 computes the next K subset of an N set. ! ! Discussion: ! ! This routine uses the revolving door method. It has no "memory". ! It simply calculates the successor of the input set, ! and will start from the beginning after the last set. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 29 March 2001 ! ! Author: ! ! FORTRAN77 original version by Albert Nijenhuis, Herbert Wilf. ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms, ! Second edition, ! Academic Press, 1978, ! ISBN 0-12-519260-6, ! LC: QA164.N54. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the size of the set from which subsets are drawn. ! N must be positive. ! ! Input, integer ( kind = 4 ) K, the size of the desired subset. K must be ! between 0 and N. ! ! Input/output, integer ( kind = 4 ) IARRAY(K). On input, the user must ! supply a subset of size K in IARRAY. That is, IARRAY must ! contain K unique numbers, in order, between 1 and N. On ! output, IARRAY(I) is the I-th element of the output subset. ! The output array is also in sorted order. ! ! Output, integer ( kind = 4 ) IN, the element of the output subset which ! was not in the input set. Each new subset differs from the ! last one by adding one element and deleting another. ! ! Output, integer ( kind = 4 ) IOUT, the element of the input subset which ! is not in the output subset. ! implicit none integer ( kind = 4 ) k integer ( kind = 4 ) iarray(k) integer ( kind = 4 ) in integer ( kind = 4 ) iout integer ( kind = 4 ) j integer ( kind = 4 ) m integer ( kind = 4 ) n if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'KSUB_NEXT2 - Fatal error!' write ( *, '(a,i8)' ) ' N = ', n write ( *, '(a)' ) ' but 0 < N is required!' stop end if if ( k < 0 .or. n < k ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'KSUB_NEXT2 - Fatal error!' write ( *, '(a,i8)' ) ' N = ', n write ( *, '(a,i8)' ) ' K = ', k write ( *, '(a)' ) ' but 0 <= K <= N is required!' stop end if j = 0 do if ( 0 < j .or. mod ( k, 2 ) == 0 ) then j = j + 1 if ( k < j ) then iarray(k) = k in = k iout = n return end if if ( iarray(j) /= j ) then iout = iarray(j) in = iout - 1 iarray(j) = in if ( j /= 1 ) then in = j - 1 iarray(j-1) = in end if return end if end if j = j + 1 m = n if ( j < k ) then m = iarray(j+1) - 1 end if if ( m /= iarray(j) ) then exit end if end do in = iarray(j) + 1 iarray(j) = in iout = in - 1 if ( j /= 1 ) then iarray(j-1) = iout iout = j - 1 end if return end subroutine legendre_set ( order, xtab, weight ) !*****************************************************************************80 ! !! LEGENDRE_SET sets abscissas and weights for Gauss-Legendre quadrature. ! ! Integration region: ! ! [ -1, 1 ] ! ! Weight function: ! ! 1.0D+00 ! ! Integral to approximate: ! ! INTEGRAL ( -1 <= X <= 1 ) F(X) dX ! ! Approximate integral: ! ! sum ( 1 <= I <= ORDER ) WEIGHT(I) * F ( XTAB(I) ) ! ! Precision: ! ! The quadrature rule will integrate exactly all polynomials up to ! X**(2*ORDER-1). ! ! Discussion: ! ! The abscissas of the rule are the zeroes of the Legendre polynomial ! P(ORDER)(X). ! ! The integral produced by a Gauss-Legendre rule is equal to the ! integral of the unique polynomial of degree ORDER-1 which ! agrees with the function at the ORDER abscissas of the rule. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 18 December 2000 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964, ! ISBN: 0-486-61272-4, ! LC: QA47.A34. ! ! Vladimir Krylov, ! Approximate Calculation of Integrals, ! Dover, 2006, ! ISBN: 0486445798, ! LC: QA311.K713. ! ! Arthur Stroud, Don Secrest, ! Gaussian Quadrature Formulas, ! Prentice Hall, 1966, ! LC: QA299.4G3S7. ! ! Daniel Zwillinger, editor, ! Standard Mathematical Tables and Formulae, ! CRC Press, 2000, ! ISBN: 0-8493-2479-3, ! LC: QA47.M315. ! ! Parameters: ! ! Input, integer ( kind = 4 ) ORDER, the order of the rule. ! ORDER must be between 1 and 20, 32 or 64. ! ! Output, real ( kind = 8 ) XTAB(ORDER), the abscissas of the rule. ! ! Output, real ( kind = 8 ) WEIGHT(ORDER), the weights of the rule. ! The weights are positive, symmetric and should sum to 2. ! implicit none integer ( kind = 4 ) order real ( kind = 8 ) xtab(order) real ( kind = 8 ) weight(order) if ( order == 1 ) then xtab(1) = 0.0D+00 weight(1) = 2.0D+00 else if ( order == 2 ) then xtab(1) = -0.577350269189625764509148780502D+00 xtab(2) = 0.577350269189625764509148780502D+00 weight(1) = 1.0D+00 weight(2) = 1.0D+00 else if ( order == 3 ) then xtab(1) = -0.774596669241483377035853079956D+00 xtab(2) = 0.0D+00 xtab(3) = 0.774596669241483377035853079956D+00 weight(1) = 5.0D+00 / 9.0D+00 weight(2) = 8.0D+00 / 9.0D+00 weight(3) = 5.0D+00 / 9.0D+00 else if ( order == 4 ) then xtab(1) = -0.861136311594052575223946488893D+00 xtab(2) = -0.339981043584856264802665759103D+00 xtab(3) = 0.339981043584856264802665759103D+00 xtab(4) = 0.861136311594052575223946488893D+00 weight(1) = 0.347854845137453857373063949222D+00 weight(2) = 0.652145154862546142626936050778D+00 weight(3) = 0.652145154862546142626936050778D+00 weight(4) = 0.347854845137453857373063949222D+00 else if ( order == 5 ) then xtab(1) = -0.906179845938663992797626878299D+00 xtab(2) = -0.538469310105683091036314420700D+00 xtab(3) = 0.0D+00 xtab(4) = 0.538469310105683091036314420700D+00 xtab(5) = 0.906179845938663992797626878299D+00 weight(1) = 0.236926885056189087514264040720D+00 weight(2) = 0.478628670499366468041291514836D+00 weight(3) = 0.568888888888888888888888888889D+00 weight(4) = 0.478628670499366468041291514836D+00 weight(5) = 0.236926885056189087514264040720D+00 else if ( order == 6 ) then xtab(1) = -0.932469514203152027812301554494D+00 xtab(2) = -0.661209386466264513661399595020D+00 xtab(3) = -0.238619186083196908630501721681D+00 xtab(4) = 0.238619186083196908630501721681D+00 xtab(5) = 0.661209386466264513661399595020D+00 xtab(6) = 0.932469514203152027812301554494D+00 weight(1) = 0.171324492379170345040296142173D+00 weight(2) = 0.360761573048138607569833513838D+00 weight(3) = 0.467913934572691047389870343990D+00 weight(4) = 0.467913934572691047389870343990D+00 weight(5) = 0.360761573048138607569833513838D+00 weight(6) = 0.171324492379170345040296142173D+00 else if ( order == 7 ) then xtab(1) = -0.949107912342758524526189684048D+00 xtab(2) = -0.741531185599394439863864773281D+00 xtab(3) = -0.405845151377397166906606412077D+00 xtab(4) = 0.0D+00 xtab(5) = 0.405845151377397166906606412077D+00 xtab(6) = 0.741531185599394439863864773281D+00 xtab(7) = 0.949107912342758524526189684048D+00 weight(1) = 0.129484966168869693270611432679D+00 weight(2) = 0.279705391489276667901467771424D+00 weight(3) = 0.381830050505118944950369775489D+00 weight(4) = 0.417959183673469387755102040816D+00 weight(5) = 0.381830050505118944950369775489D+00 weight(6) = 0.279705391489276667901467771424D+00 weight(7) = 0.129484966168869693270611432679D+00 else if ( order == 8 ) then xtab(1) = -0.960289856497536231683560868569D+00 xtab(2) = -0.796666477413626739591553936476D+00 xtab(3) = -0.525532409916328985817739049189D+00 xtab(4) = -0.183434642495649804939476142360D+00 xtab(5) = 0.183434642495649804939476142360D+00 xtab(6) = 0.525532409916328985817739049189D+00 xtab(7) = 0.796666477413626739591553936476D+00 xtab(8) = 0.960289856497536231683560868569D+00 weight(1) = 0.101228536290376259152531354310D+00 weight(2) = 0.222381034453374470544355994426D+00 weight(3) = 0.313706645877887287337962201987D+00 weight(4) = 0.362683783378361982965150449277D+00 weight(5) = 0.362683783378361982965150449277D+00 weight(6) = 0.313706645877887287337962201987D+00 weight(7) = 0.222381034453374470544355994426D+00 weight(8) = 0.101228536290376259152531354310D+00 else if ( order == 9 ) then xtab(1) = -0.968160239507626089835576202904D+00 xtab(2) = -0.836031107326635794299429788070D+00 xtab(3) = -0.613371432700590397308702039341D+00 xtab(4) = -0.324253423403808929038538014643D+00 xtab(5) = 0.0D+00 xtab(6) = 0.324253423403808929038538014643D+00 xtab(7) = 0.613371432700590397308702039341D+00 xtab(8) = 0.836031107326635794299429788070D+00 xtab(9) = 0.968160239507626089835576202904D+00 weight(1) = 0.812743883615744119718921581105D-01 weight(2) = 0.180648160694857404058472031243D+00 weight(3) = 0.260610696402935462318742869419D+00 weight(4) = 0.312347077040002840068630406584D+00 weight(5) = 0.330239355001259763164525069287D+00 weight(6) = 0.312347077040002840068630406584D+00 weight(7) = 0.260610696402935462318742869419D+00 weight(8) = 0.180648160694857404058472031243D+00 weight(9) = 0.812743883615744119718921581105D-01 else if ( order == 10 ) then xtab(1) = -0.973906528517171720077964012084D+00 xtab(2) = -0.865063366688984510732096688423D+00 xtab(3) = -0.679409568299024406234327365115D+00 xtab(4) = -0.433395394129247190799265943166D+00 xtab(5) = -0.148874338981631210884826001130D+00 xtab(6) = 0.148874338981631210884826001130D+00 xtab(7) = 0.433395394129247190799265943166D+00 xtab(8) = 0.679409568299024406234327365115D+00 xtab(9) = 0.865063366688984510732096688423D+00 xtab(10) = 0.973906528517171720077964012084D+00 weight(1) = 0.666713443086881375935688098933D-01 weight(2) = 0.149451349150580593145776339658D+00 weight(3) = 0.219086362515982043995534934228D+00 weight(4) = 0.269266719309996355091226921569D+00 weight(5) = 0.295524224714752870173892994651D+00 weight(6) = 0.295524224714752870173892994651D+00 weight(7) = 0.269266719309996355091226921569D+00 weight(8) = 0.219086362515982043995534934228D+00 weight(9) = 0.149451349150580593145776339658D+00 weight(10) = 0.666713443086881375935688098933D-01 else if ( order == 11 ) then xtab(1) = -0.978228658146056992803938001123D+00 xtab(2) = -0.887062599768095299075157769304D+00 xtab(3) = -0.730152005574049324093416252031D+00 xtab(4) = -0.519096129206811815925725669459D+00 xtab(5) = -0.269543155952344972331531985401D+00 xtab(6) = 0.0D+00 xtab(7) = 0.269543155952344972331531985401D+00 xtab(8) = 0.519096129206811815925725669459D+00 xtab(9) = 0.730152005574049324093416252031D+00 xtab(10) = 0.887062599768095299075157769304D+00 xtab(11) = 0.978228658146056992803938001123D+00 weight(1) = 0.556685671161736664827537204425D-01 weight(2) = 0.125580369464904624634694299224D+00 weight(3) = 0.186290210927734251426097641432D+00 weight(4) = 0.233193764591990479918523704843D+00 weight(5) = 0.262804544510246662180688869891D+00 weight(6) = 0.272925086777900630714483528336D+00 weight(7) = 0.262804544510246662180688869891D+00 weight(8) = 0.233193764591990479918523704843D+00 weight(9) = 0.186290210927734251426097641432D+00 weight(10) = 0.125580369464904624634694299224D+00 weight(11) = 0.556685671161736664827537204425D-01 else if ( order == 12 ) then xtab(1) = -0.981560634246719250690549090149D+00 xtab(2) = -0.904117256370474856678465866119D+00 xtab(3) = -0.769902674194304687036893833213D+00 xtab(4) = -0.587317954286617447296702418941D+00 xtab(5) = -0.367831498998180193752691536644D+00 xtab(6) = -0.125233408511468915472441369464D+00 xtab(7) = 0.125233408511468915472441369464D+00 xtab(8) = 0.367831498998180193752691536644D+00 xtab(9) = 0.587317954286617447296702418941D+00 xtab(10) = 0.769902674194304687036893833213D+00 xtab(11) = 0.904117256370474856678465866119D+00 xtab(12) = 0.981560634246719250690549090149D+00 weight(1) = 0.471753363865118271946159614850D-01 weight(2) = 0.106939325995318430960254718194D+00 weight(3) = 0.160078328543346226334652529543D+00 weight(4) = 0.203167426723065921749064455810D+00 weight(5) = 0.233492536538354808760849898925D+00 weight(6) = 0.249147045813402785000562436043D+00 weight(7) = 0.249147045813402785000562436043D+00 weight(8) = 0.233492536538354808760849898925D+00 weight(9) = 0.203167426723065921749064455810D+00 weight(10) = 0.160078328543346226334652529543D+00 weight(11) = 0.106939325995318430960254718194D+00 weight(12) = 0.471753363865118271946159614850D-01 else if ( order == 13 ) then xtab(1) = -0.984183054718588149472829448807D+00 xtab(2) = -0.917598399222977965206547836501D+00 xtab(3) = -0.801578090733309912794206489583D+00 xtab(4) = -0.642349339440340220643984606996D+00 xtab(5) = -0.448492751036446852877912852128D+00 xtab(6) = -0.230458315955134794065528121098D+00 xtab(7) = 0.0D+00 xtab(8) = 0.230458315955134794065528121098D+00 xtab(9) = 0.448492751036446852877912852128D+00 xtab(10) = 0.642349339440340220643984606996D+00 xtab(11) = 0.801578090733309912794206489583D+00 xtab(12) = 0.917598399222977965206547836501D+00 xtab(13) = 0.984183054718588149472829448807D+00 weight(1) = 0.404840047653158795200215922010D-01 weight(2) = 0.921214998377284479144217759538D-01 weight(3) = 0.138873510219787238463601776869D+00 weight(4) = 0.178145980761945738280046691996D+00 weight(5) = 0.207816047536888502312523219306D+00 weight(6) = 0.226283180262897238412090186040D+00 weight(7) = 0.232551553230873910194589515269D+00 weight(8) = 0.226283180262897238412090186040D+00 weight(9) = 0.207816047536888502312523219306D+00 weight(10) = 0.178145980761945738280046691996D+00 weight(11) = 0.138873510219787238463601776869D+00 weight(12) = 0.921214998377284479144217759538D-01 weight(13) = 0.404840047653158795200215922010D-01 else if ( order == 14 ) then xtab(1) = -0.986283808696812338841597266704D+00 xtab(2) = -0.928434883663573517336391139378D+00 xtab(3) = -0.827201315069764993189794742650D+00 xtab(4) = -0.687292904811685470148019803019D+00 xtab(5) = -0.515248636358154091965290718551D+00 xtab(6) = -0.319112368927889760435671824168D+00 xtab(7) = -0.108054948707343662066244650220D+00 xtab(8) = 0.108054948707343662066244650220D+00 xtab(9) = 0.319112368927889760435671824168D+00 xtab(10) = 0.515248636358154091965290718551D+00 xtab(11) = 0.687292904811685470148019803019D+00 xtab(12) = 0.827201315069764993189794742650D+00 xtab(13) = 0.928434883663573517336391139378D+00 xtab(14) = 0.986283808696812338841597266704D+00 weight(1) = 0.351194603317518630318328761382D-01 weight(2) = 0.801580871597602098056332770629D-01 weight(3) = 0.121518570687903184689414809072D+00 weight(4) = 0.157203167158193534569601938624D+00 weight(5) = 0.185538397477937813741716590125D+00 weight(6) = 0.205198463721295603965924065661D+00 weight(7) = 0.215263853463157790195876443316D+00 weight(8) = 0.215263853463157790195876443316D+00 weight(9) = 0.205198463721295603965924065661D+00 weight(10) = 0.185538397477937813741716590125D+00 weight(11) = 0.157203167158193534569601938624D+00 weight(12) = 0.121518570687903184689414809072D+00 weight(13) = 0.801580871597602098056332770629D-01 weight(14) = 0.351194603317518630318328761382D-01 else if ( order == 15 ) then xtab(1) = -0.987992518020485428489565718587D+00 xtab(2) = -0.937273392400705904307758947710D+00 xtab(3) = -0.848206583410427216200648320774D+00 xtab(4) = -0.724417731360170047416186054614D+00 xtab(5) = -0.570972172608538847537226737254D+00 xtab(6) = -0.394151347077563369897207370981D+00 xtab(7) = -0.201194093997434522300628303395D+00 xtab(8) = 0.0D+00 xtab(9) = 0.201194093997434522300628303395D+00 xtab(10) = 0.394151347077563369897207370981D+00 xtab(11) = 0.570972172608538847537226737254D+00 xtab(12) = 0.724417731360170047416186054614D+00 xtab(13) = 0.848206583410427216200648320774D+00 xtab(14) = 0.937273392400705904307758947710D+00 xtab(15) = 0.987992518020485428489565718587D+00 weight(1) = 0.307532419961172683546283935772D-01 weight(2) = 0.703660474881081247092674164507D-01 weight(3) = 0.107159220467171935011869546686D+00 weight(4) = 0.139570677926154314447804794511D+00 weight(5) = 0.166269205816993933553200860481D+00 weight(6) = 0.186161000015562211026800561866D+00 weight(7) = 0.198431485327111576456118326444D+00 weight(8) = 0.202578241925561272880620199968D+00 weight(9) = 0.198431485327111576456118326444D+00 weight(10) = 0.186161000015562211026800561866D+00 weight(11) = 0.166269205816993933553200860481D+00 weight(12) = 0.139570677926154314447804794511D+00 weight(13) = 0.107159220467171935011869546686D+00 weight(14) = 0.703660474881081247092674164507D-01 weight(15) = 0.307532419961172683546283935772D-01 else if ( order == 16 ) then xtab(1) = -0.989400934991649932596154173450D+00 xtab(2) = -0.944575023073232576077988415535D+00 xtab(3) = -0.865631202387831743880467897712D+00 xtab(4) = -0.755404408355003033895101194847D+00 xtab(5) = -0.617876244402643748446671764049D+00 xtab(6) = -0.458016777657227386342419442984D+00 xtab(7) = -0.281603550779258913230460501460D+00 xtab(8) = -0.950125098376374401853193354250D-01 xtab(9) = 0.950125098376374401853193354250D-01 xtab(10) = 0.281603550779258913230460501460D+00 xtab(11) = 0.458016777657227386342419442984D+00 xtab(12) = 0.617876244402643748446671764049D+00 xtab(13) = 0.755404408355003033895101194847D+00 xtab(14) = 0.865631202387831743880467897712D+00 xtab(15) = 0.944575023073232576077988415535D+00 xtab(16) = 0.989400934991649932596154173450D+00 weight(1) = 0.271524594117540948517805724560D-01 weight(2) = 0.622535239386478928628438369944D-01 weight(3) = 0.951585116824927848099251076022D-01 weight(4) = 0.124628971255533872052476282192D+00 weight(5) = 0.149595988816576732081501730547D+00 weight(6) = 0.169156519395002538189312079030D+00 weight(7) = 0.182603415044923588866763667969D+00 weight(8) = 0.189450610455068496285396723208D+00 weight(9) = 0.189450610455068496285396723208D+00 weight(10) = 0.182603415044923588866763667969D+00 weight(11) = 0.169156519395002538189312079030D+00 weight(12) = 0.149595988816576732081501730547D+00 weight(13) = 0.124628971255533872052476282192D+00 weight(14) = 0.951585116824927848099251076022D-01 weight(15) = 0.622535239386478928628438369944D-01 weight(16) = 0.271524594117540948517805724560D-01 else if ( order == 17 ) then xtab(1) = -0.990575475314417335675434019941D+00 xtab(2) = -0.950675521768767761222716957896D+00 xtab(3) = -0.880239153726985902122955694488D+00 xtab(4) = -0.781514003896801406925230055520D+00 xtab(5) = -0.657671159216690765850302216643D+00 xtab(6) = -0.512690537086476967886246568630D+00 xtab(7) = -0.351231763453876315297185517095D+00 xtab(8) = -0.178484181495847855850677493654D+00 xtab(9) = 0.0D+00 xtab(10) = 0.178484181495847855850677493654D+00 xtab(11) = 0.351231763453876315297185517095D+00 xtab(12) = 0.512690537086476967886246568630D+00 xtab(13) = 0.657671159216690765850302216643D+00 xtab(14) = 0.781514003896801406925230055520D+00 xtab(15) = 0.880239153726985902122955694488D+00 xtab(16) = 0.950675521768767761222716957896D+00 xtab(17) = 0.990575475314417335675434019941D+00 weight(1) = 0.241483028685479319601100262876D-01 weight(2) = 0.554595293739872011294401653582D-01 weight(3) = 0.850361483171791808835353701911D-01 weight(4) = 0.111883847193403971094788385626D+00 weight(5) = 0.135136368468525473286319981702D+00 weight(6) = 0.154045761076810288081431594802D+00 weight(7) = 0.168004102156450044509970663788D+00 weight(8) = 0.176562705366992646325270990113D+00 weight(9) = 0.179446470356206525458265644262D+00 weight(10) = 0.176562705366992646325270990113D+00 weight(11) = 0.168004102156450044509970663788D+00 weight(12) = 0.154045761076810288081431594802D+00 weight(13) = 0.135136368468525473286319981702D+00 weight(14) = 0.111883847193403971094788385626D+00 weight(15) = 0.850361483171791808835353701911D-01 weight(16) = 0.554595293739872011294401653582D-01 weight(17) = 0.241483028685479319601100262876D-01 else if ( order == 18 ) then xtab(1) = -0.991565168420930946730016004706D+00 xtab(2) = -0.955823949571397755181195892930D+00 xtab(3) = -0.892602466497555739206060591127D+00 xtab(4) = -0.803704958972523115682417455015D+00 xtab(5) = -0.691687043060353207874891081289D+00 xtab(6) = -0.559770831073947534607871548525D+00 xtab(7) = -0.411751161462842646035931793833D+00 xtab(8) = -0.251886225691505509588972854878D+00 xtab(9) = -0.847750130417353012422618529358D-01 xtab(10) = 0.847750130417353012422618529358D-01 xtab(11) = 0.251886225691505509588972854878D+00 xtab(12) = 0.411751161462842646035931793833D+00 xtab(13) = 0.559770831073947534607871548525D+00 xtab(14) = 0.691687043060353207874891081289D+00 xtab(15) = 0.803704958972523115682417455015D+00 xtab(16) = 0.892602466497555739206060591127D+00 xtab(17) = 0.955823949571397755181195892930D+00 xtab(18) = 0.991565168420930946730016004706D+00 weight(1) = 0.216160135264833103133427102665D-01 weight(2) = 0.497145488949697964533349462026D-01 weight(3) = 0.764257302548890565291296776166D-01 weight(4) = 0.100942044106287165562813984925D+00 weight(5) = 0.122555206711478460184519126800D+00 weight(6) = 0.140642914670650651204731303752D+00 weight(7) = 0.154684675126265244925418003836D+00 weight(8) = 0.164276483745832722986053776466D+00 weight(9) = 0.169142382963143591840656470135D+00 weight(10) = 0.169142382963143591840656470135D+00 weight(11) = 0.164276483745832722986053776466D+00 weight(12) = 0.154684675126265244925418003836D+00 weight(13) = 0.140642914670650651204731303752D+00 weight(14) = 0.122555206711478460184519126800D+00 weight(15) = 0.100942044106287165562813984925D+00 weight(16) = 0.764257302548890565291296776166D-01 weight(17) = 0.497145488949697964533349462026D-01 weight(18) = 0.216160135264833103133427102665D-01 else if ( order == 19 ) then xtab(1) = -0.992406843843584403189017670253D+00 xtab(2) = -0.960208152134830030852778840688D+00 xtab(3) = -0.903155903614817901642660928532D+00 xtab(4) = -0.822714656537142824978922486713D+00 xtab(5) = -0.720966177335229378617095860824D+00 xtab(6) = -0.600545304661681023469638164946D+00 xtab(7) = -0.464570741375960945717267148104D+00 xtab(8) = -0.316564099963629831990117328850D+00 xtab(9) = -0.160358645640225375868096115741D+00 xtab(10) = 0.0D+00 xtab(11) = 0.160358645640225375868096115741D+00 xtab(12) = 0.316564099963629831990117328850D+00 xtab(13) = 0.464570741375960945717267148104D+00 xtab(14) = 0.600545304661681023469638164946D+00 xtab(15) = 0.720966177335229378617095860824D+00 xtab(16) = 0.822714656537142824978922486713D+00 xtab(17) = 0.903155903614817901642660928532D+00 xtab(18) = 0.960208152134830030852778840688D+00 xtab(19) = 0.992406843843584403189017670253D+00 weight(1) = 0.194617882297264770363120414644D-01 weight(2) = 0.448142267656996003328381574020D-01 weight(3) = 0.690445427376412265807082580060D-01 weight(4) = 0.914900216224499994644620941238D-01 weight(5) = 0.111566645547333994716023901682D+00 weight(6) = 0.128753962539336227675515784857D+00 weight(7) = 0.142606702173606611775746109442D+00 weight(8) = 0.152766042065859666778855400898D+00 weight(9) = 0.158968843393954347649956439465D+00 weight(10) = 0.161054449848783695979163625321D+00 weight(11) = 0.158968843393954347649956439465D+00 weight(12) = 0.152766042065859666778855400898D+00 weight(13) = 0.142606702173606611775746109442D+00 weight(14) = 0.128753962539336227675515784857D+00 weight(15) = 0.111566645547333994716023901682D+00 weight(16) = 0.914900216224499994644620941238D-01 weight(17) = 0.690445427376412265807082580060D-01 weight(18) = 0.448142267656996003328381574020D-01 weight(19) = 0.194617882297264770363120414644D-01 else if ( order == 20 ) then xtab(1) = -0.993128599185094924786122388471D+00 xtab(2) = -0.963971927277913791267666131197D+00 xtab(3) = -0.912234428251325905867752441203D+00 xtab(4) = -0.839116971822218823394529061702D+00 xtab(5) = -0.746331906460150792614305070356D+00 xtab(6) = -0.636053680726515025452836696226D+00 xtab(7) = -0.510867001950827098004364050955D+00 xtab(8) = -0.373706088715419560672548177025D+00 xtab(9) = -0.227785851141645078080496195369D+00 xtab(10) = -0.765265211334973337546404093988D-01 xtab(11) = 0.765265211334973337546404093988D-01 xtab(12) = 0.227785851141645078080496195369D+00 xtab(13) = 0.373706088715419560672548177025D+00 xtab(14) = 0.510867001950827098004364050955D+00 xtab(15) = 0.636053680726515025452836696226D+00 xtab(16) = 0.746331906460150792614305070356D+00 xtab(17) = 0.839116971822218823394529061702D+00 xtab(18) = 0.912234428251325905867752441203D+00 xtab(19) = 0.963971927277913791267666131197D+00 xtab(20) = 0.993128599185094924786122388471D+00 weight(1) = 0.176140071391521183118619623519D-01 weight(2) = 0.406014298003869413310399522749D-01 weight(3) = 0.626720483341090635695065351870D-01 weight(4) = 0.832767415767047487247581432220D-01 weight(5) = 0.101930119817240435036750135480D+00 weight(6) = 0.118194531961518417312377377711D+00 weight(7) = 0.131688638449176626898494499748D+00 weight(8) = 0.142096109318382051329298325067D+00 weight(9) = 0.149172986472603746787828737002D+00 weight(10) = 0.152753387130725850698084331955D+00 weight(11) = 0.152753387130725850698084331955D+00 weight(12) = 0.149172986472603746787828737002D+00 weight(13) = 0.142096109318382051329298325067D+00 weight(14) = 0.131688638449176626898494499748D+00 weight(15) = 0.118194531961518417312377377711D+00 weight(16) = 0.101930119817240435036750135480D+00 weight(17) = 0.832767415767047487247581432220D-01 weight(18) = 0.626720483341090635695065351870D-01 weight(19) = 0.406014298003869413310399522749D-01 weight(20) = 0.176140071391521183118619623519D-01 else if ( order == 32 ) then xtab(1) = -0.997263861849481563544981128665D+00 xtab(2) = -0.985611511545268335400175044631D+00 xtab(3) = -0.964762255587506430773811928118D+00 xtab(4) = -0.934906075937739689170919134835D+00 xtab(5) = -0.896321155766052123965307243719D+00 xtab(6) = -0.849367613732569970133693004968D+00 xtab(7) = -0.794483795967942406963097298970D+00 xtab(8) = -0.732182118740289680387426665091D+00 xtab(9) = -0.663044266930215200975115168663D+00 xtab(10) = -0.587715757240762329040745476402D+00 xtab(11) = -0.506899908932229390023747474378D+00 xtab(12) = -0.421351276130635345364119436172D+00 xtab(13) = -0.331868602282127649779916805730D+00 xtab(14) = -0.239287362252137074544603209166D+00 xtab(15) = -0.144471961582796493485186373599D+00 xtab(16) = -0.483076656877383162348125704405D-01 xtab(17) = 0.483076656877383162348125704405D-01 xtab(18) = 0.144471961582796493485186373599D+00 xtab(19) = 0.239287362252137074544603209166D+00 xtab(20) = 0.331868602282127649779916805730D+00 xtab(21) = 0.421351276130635345364119436172D+00 xtab(22) = 0.506899908932229390023747474378D+00 xtab(23) = 0.587715757240762329040745476402D+00 xtab(24) = 0.663044266930215200975115168663D+00 xtab(25) = 0.732182118740289680387426665091D+00 xtab(26) = 0.794483795967942406963097298970D+00 xtab(27) = 0.849367613732569970133693004968D+00 xtab(28) = 0.896321155766052123965307243719D+00 xtab(29) = 0.934906075937739689170919134835D+00 xtab(30) = 0.964762255587506430773811928118D+00 xtab(31) = 0.985611511545268335400175044631D+00 xtab(32) = 0.997263861849481563544981128665D+00 weight(1) = 0.70186100094700966004