program main

!*****************************************************************************80
!
!! MAIN is the main program for STROUD_PRB.
!
!  Discussion:
!
!    STROUD_PRB tests the routines in the STROUD library.
!
!  Modified:
!
!    24 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'STROUD_PRB'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Tests for the routines in the STROUD library,'
  write ( *, '(a)' ) '  which approximate integrals over N-dimensional regions.'

  call test01
  call test02
  call test03
  call test04
  call test045
  call test05
  call test052
  call test054
  call test07
  call test08
  call test085
  call test09

  call test10
  call test11
  call test12
  call test13
  call test14
  call test15
  call test16
  call test17
  call test18
  call test19

  call test20
  call test205
  call test21
  call test215
  call test22
  call test23
  call test24
  call test25
  call test255
  call test26
  call test27
  call test28
  call test29

  call test30
  call test31
  call test32
  call test322
  call test324
  call test326
  call test33
  call test335
  call test34
  call test345
  call test35
  call test36
  call test37
  call test38
  call test39

  call test40
  call test41
  call test42
  call test425
  call test43
  call test44
  call test45
  call test46
  call test47
  call test48
  call test49

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'STROUD_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine test01

!*****************************************************************************80
!
!! TEST01 tests BALL_F1_ND, BALL_F3_ND.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: n_max = 3

  real      ( kind = 8 ) ball_volume_nd
  real      ( kind = 8 ), dimension ( n_max ) :: center = (/ &
    1.0D+00, -1.0D+00, 2.0D+00 /)
  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) function_nd_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) r

  r = 2.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  For integrals in a ball in ND:'
  write ( *, '(a)' ) '  BALL_F1_ND approximates the integral;'
  write ( *, '(a)' ) '  BALL_F3_ND approximates the integral.'
  write ( *, '(a)' ) ' '

  do n = 2, n_max

    do i = 1, n
      center(i) = real ( i, kind = 8 )
    end do

    write ( *, '(a)' ) ' '
    write ( *, '(a,i4)' ) '  Spatial dimension N = ', n
    write ( *, '(a)' ) '  Ball center:'
    write ( *, '(2x,3f10.4)' ) center(1:n)
    write ( *, '(a,g14.6)') '  Ball radius = ', r
    write ( *, '(a,g14.6)' ) '  Ball volume = ', ball_volume_nd ( n, r )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '    Rule:      F1          F3'
    write ( *, '(a)' ) '    F(X)'
    write ( *, '(a)' ) ' '

    num = function_nd_num ( )

    do i = 1, num

      call function_nd_set ( 'SET', i )

      call ball_f1_nd ( function_nd, n, center, r, result1 )
      call ball_f3_nd ( function_nd, n, center, r, result2 )
      write ( *, '(2x,a7,2f14.8)' ) function_nd_name(i), result1, result2

    end do

  end do

  return
end
subroutine test02

!*****************************************************************************80
!
!! TEST02 tests BALL_MONOMIAL_ND.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 3
  integer   ( kind = 4 ), parameter :: test_num = 4

  real      ( kind = 8 ) ball_monomial_nd
  real      ( kind = 8 ) ball_volume_nd
  real      ( kind = 8 ), dimension ( dim_num ) :: center = (/ &
    0.0D+00, 0.0D+00, 0.0D+00 /)
  real      ( kind = 8 ), external :: mono_000_3d
  real      ( kind = 8 ), external :: mono_111_3d
  real      ( kind = 8 ), external :: mono_202_3d
  real      ( kind = 8 ), external :: mono_422_3d
  integer   ( kind = 4 ) p(dim_num)
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) :: r = 2.0D+00
  character ( len = 10 ) string
  integer   ( kind = 4 ) test

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  For the integral of a monomial in a ball in ND:'
  write ( *, '(a)' ) '  BALL_MONOMIAL_ND approximates the integral.'
  write ( *, '(a)' ) '  BALL_F1_ND, which can handle general integrands,'
  write ( *, '(a)' ) '    will be used for comparison.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Spatial dimension N = ', dim_num
  write ( *, '(a,g14.6)' ) '  Ball radius = ', r
  write ( *, '(a,g14.6)' ) '  Ball volume = ', &
    ball_volume_nd ( dim_num, r )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    Rule:     MONOMIAL    F1'
  write ( *, '(a)' ) '    F(X)'
  write ( *, '(a)' ) ' '

  do test = 1, test_num

    if ( test == 1 ) then
      string = '1'
      p = (/ 0, 0, 0 /)
      call ball_f1_nd ( mono_000_3d, dim_num, center, r, result2 )
    else if ( test == 2 ) then
      string = 'xyz'
      p = (/ 1, 1, 1 /)
      call ball_f1_nd ( mono_111_3d, dim_num, center, r, result2 )
    else if ( test == 3 ) then
      string = 'x^2z^2'
      p = (/ 2, 0, 2 /)
      call ball_f1_nd ( mono_202_3d, dim_num, center, r, result2 )
    else if ( test == 4 ) then
      string = 'x^4y^2z^2'
      p = (/ 4, 2, 2 /)
      call ball_f1_nd ( mono_422_3d, dim_num, center, r, result2 )
    end if

    result1 = ball_monomial_nd ( dim_num, p, r )

    write ( *, '(2x,a10,2f14.8)' ) string, result1, result2

  end do

  return
end
subroutine test03

!*****************************************************************************80
!
!! TEST03 tests BALL_UNIT_**_3D.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ) ball_unit_volume_nd
  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) result3

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  For integrals in the unit ball in 3D:'
  write ( *, '(a)' ) '  BALL_UNIT_07_3D uses a formula of degree 7;'
  write ( *, '(a)' ) '  BALL_UNIT_14_3D uses a formula of degree 14;'
  write ( *, '(a)' ) '  BALL_UNIT_15_3D uses a formula of degree 15.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Unit ball volume = ', ball_unit_volume_nd ( 3 )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    Rule:      #7             #14           #15'
  write ( *, '(a)' ) '    F(X)'
  write ( *, '(a)' ) ' '
 
  num = function_3d_num ( )

  do i = 1, num

    call function_3d_set ( 'SET', i )

    call ball_unit_07_3d ( function_3d, result1 )
    call ball_unit_14_3d ( function_3d, result2 )
    call ball_unit_15_3d ( function_3d, result3 )

    write ( *, '(2x,a7,2x,g14.6,2x,g14.6,2x,g14.6)' ) &
      function_3d_name(i), result1, result2, result3

  end do
 
  return
end
subroutine test04

!*****************************************************************************80
!
!! TEST04 tests BALL_UNIT_F1_ND, BALL_UNIT_F3_ND.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: n_max = 3

  real      ( kind = 8 ) ball_unit_volume_nd
  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) function_nd_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  For integrals inside the unit ball in ND:'
  write ( *, '(a)' ) '  BALL_UNIT_F1_ND approximates the integral;'
  write ( *, '(a)' ) '  BALL_UNIT_F3_ND approximates the integral.'
  write ( *, '(a)' ) ' '
 
  do n = 2, n_max
 
    write ( *, '(a)' ) ' '
    write ( *, '(a,i4)' ) '  Spatial dimension N = ', n
    write ( *, '(a,g14.6)' ) '  Unit ball volume = ', ball_unit_volume_nd ( n )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '    Rule:      F1          F3'
    write ( *, '(a)' ) '    F(X)'
    write ( *, '(a)' ) ' '
 
    num = function_nd_num ( )

    do i = 1, num

      call function_nd_set ( 'SET', i )

      call ball_unit_f1_nd ( function_nd, n, result1 )
      call ball_unit_f3_nd ( function_nd, n, result2 )
      write ( *, '(2x,a7,2f14.8)' ) function_nd_name(i), result1, result2
 
    end do
 
  end do
 
  return
end
subroutine test045

!*****************************************************************************80
!
!! TEST045 tests BALL_UNIT_VOLUME_3D.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ) ball_unit_volume_3d
  real      ( kind = 8 ) ball_unit_volume_nd
  integer   ( kind = 4 ), parameter :: dim_num = 3

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST045'
  write ( *, '(a)' ) '  In 3 dimensions:'
  write ( *, '(a)' ) '  BALL_UNIT_VOLUME_3D gets the volume of the unit ball.'
  write ( *, '(a)' ) '  BALL_UNIT_VOLUME_ND will be called for comparison.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    N    Volume    Method'
  write ( *, '(a)' ) ' '

  write ( *, '(2x,i3,f14.8,2x,a)' ) dim_num, &
    ball_unit_volume_3d (   ), 'BALL_UNIT_VOLUME_3D'
  write ( *, '(2x,i3,f14.8,2x,a)' ) dim_num, &
    ball_unit_volume_nd ( dim_num ), 'BALL_UNIT_VOLUME_ND'

  return
end
subroutine test05

!*****************************************************************************80
!
!! TEST05 tests BALL_UNIT_VOLUME_ND.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none
 
  real      ( kind = 8 ) ball_unit_volume_nd
  integer   ( kind = 4 ) dim_num

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  BALL_UNIT_VOLUME_ND computes the volume '
  write ( *, '(a)' ) '    of the unit ball in ND.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    N    Volume'
  write ( *, '(a)' ) ' '

  do dim_num = 2, 10
    write ( *, '(2x,i3,f14.8)' ) dim_num, &
      ball_unit_volume_nd ( dim_num )
  end do

  return
end
subroutine test052

!*****************************************************************************80
!
!! TEST052 tests BALL_VOLUME_3D.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ) ball_volume_3d
  real      ( kind = 8 ) ball_volume_nd
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ), parameter :: n = 3
  real      ( kind = 8 ) r

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST052'
  write ( *, '(a)' ) '  In 3 dimensions:'
  write ( *, '(a)' ) '  BALL_VOLUME_3D computes the volume of a unit ball.'
  write ( *, '(a)' ) '  BALL_VOLUME_ND will be called for comparison.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    N    R      Volume    Method'
  write ( *, '(a)' ) ' '

  r = 1.0D+00

  do i = 1, 3

    write ( *, '(2x,i3,2f14.8,2x,a)' ) &
      n, r, ball_volume_3d (    r ), 'BALL_VOLUME_3D'

    write ( *, '(2x,i3,2f14.8,2x,a)' ) &
      n, r, ball_volume_nd ( n, r ), 'BALL_VOLUME_ND'

    r = r * 2.0D+00

  end do

  return
end
subroutine test054

!*****************************************************************************80
!
!! TEST054 tests BALL_VOLUME_ND.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ) ball_volume_nd
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) n
  real      ( kind = 8 ) r

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST054'
  write ( *, '(a)' ) '  BALL_UNIT_VOLUME_ND computes the volume of '
  write ( *, '(a)' ) '    the unit ball in N dimensions.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    N        R      Volume'
  write ( *, '(a)' ) ' '

  do n = 2, 10
    r = 0.5D+00
    do i = 1, 3
      write ( *, '(2x,i3,2f14.8)' ) n, r, ball_volume_nd ( n, r )
      r = r * 2.0D+00
    end do
  end do

  return
end
subroutine test07

!*****************************************************************************80
!
!! TEST07 tests CIRCLE_ANNULUS.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2
  integer   ( kind = 4 ), parameter :: test_num = 2

  real      ( kind = 8 ) area
  real      ( kind = 8 ) center(dim_num)
  real      ( kind = 8 ), dimension ( dim_num, test_num ) :: center_test = reshape ( (/ &
    0.0D+00, 0.0D+00, &
    0.0D+00, 0.0D+00 /), (/ dim_num, test_num /) )
  real      ( kind = 8 ) circle_annulus_area_2d
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) nr
  real      ( kind = 8 ) radius1
  real      ( kind = 8 ), dimension ( test_num ) :: radius1_test = (/ &
    0.0D+00, 1.0D+00 /)
  real      ( kind = 8 ) radius2
  real      ( kind = 8 ), dimension ( test_num ) :: radius2_test = (/ &
    1.0D+00, 2.0D+00 /)
  real      ( kind = 8 ) result

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  CIRCLE_ANNULUS estimates integrals '
  write ( *, '(a)' ) '    in a circular annulus.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) &
    '        F       CENTER         Radius1   Radius2   NR  Result'
  write ( *, '(a)' ) ' '

  do i = 1, test_num

    center(1:dim_num) = center_test(1:dim_num,i)

    radius1 = radius1_test(i)
    radius2 = radius2_test(i)

    area = circle_annulus_area_2d ( radius1, radius2 )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a7,4f10.6,2x,f10.6)' ) &
      '   Area', center(1:dim_num), radius1, radius2, area

    num = function_2d_num ( )

    do j = 1, num

      call function_2d_set ( 'SET', j )

      do nr = 1, 4
        call circle_annulus ( function_2d, center, radius1, radius2, nr, result )
        write ( *, '(2x,a7,4f10.6,i2,f10.6)' ) &
          function_2d_name(j), center(1:dim_num), radius1, radius2, nr, result
      end do

    end do

  end do

  return
end
subroutine test08

!*****************************************************************************80
!
!! TEST08 tests CIRCLE_ANNULUS, CIRCLE_RT_SET, CIRCLE_RT_SUM.
!
!  Modified:
!
!    07 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2
  integer   ( kind = 4 ), parameter :: test_num = 3

  real      ( kind = 8 ) area
  real      ( kind = 8 ) center(dim_num)
  real      ( kind = 8 ), dimension ( dim_num, test_num ) :: center_test = reshape ( (/ &
    0.0D+00, 0.0D+00, &
    0.0D+00, 0.0D+00, &
    0.0D+00, 0.0D+00 /), (/ dim_num, test_num /) )
  real      ( kind = 8 ) circle_annulus_area_2d
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) nc
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) nr
  integer   ( kind = 4 ) nr2
  integer   ( kind = 4 ) nt
  real      ( kind = 8 ) ra(5)
  real      ( kind = 8 ) radius1
  real      ( kind = 8 ), dimension ( test_num ) :: radius1_test = (/ &
    0.0D+00, 1.0D+00, 1.0D+00 /)
  real      ( kind = 8 ) radius2
  real      ( kind = 8 ), dimension ( test_num ) :: radius2_test = (/ &
    1.0D+00, 2.0D+00, 3.0D+00 /)
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) result3
  real      ( kind = 8 ) rw(5)
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ) ta(20)
  real      ( kind = 8 ) tw(20)
  real      ( kind = 8 ) zw

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  CIRCLE_ANNULUS estimates integrals in a '
  write ( *, '(a)' ) '    circular annulus.'
  write ( *, '(a)' ) '  CIRCLE_RT_SET sets up a rule for a circle;'
  write ( *, '(a)' ) '  CIRCLE_RT_SUM applies the rule.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  RESULT1 = CIRCLE_ANNULUS result.'
  write ( *, '(a)' ) '  RESULT2 = Difference of two CIRCLE_RT_SUM results.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '        F      CENTER       Radius1   Radius2   Result1 Result2'
  write ( *, '(a)' ) ' '

  do i = 1, test_num

    center(1:dim_num) = center_test(1:dim_num,i)

    radius1 = radius1_test(i)
    radius2 = radius2_test(i)

    area = circle_annulus_area_2d ( radius1, radius2 )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a7,5f11.5)' ) &
      '   Area', center(1:dim_num), radius1, radius2, area

    rule = 9
    call circle_rt_size ( rule, nr2, nt, nc )
    call circle_rt_set ( rule, nr2, nt, nc, ra, rw, ta, tw, zw )

    num = function_2d_num ( )

    do j = 1, num

      call function_2d_set ( 'SET', j )

      nr = 5
      call circle_annulus ( function_2d, center, radius1, radius2, nr, result1 )

      call circle_rt_sum ( function_2d, center, radius1, nr2, ra, rw, nt, ta, tw, &
        zw, result2 )

      call circle_rt_sum ( function_2d, center, radius2, nr2, ra, rw, nt, ta, tw, &
        zw, result3 )

      write ( *, '(2x,a7,6f11.5)' ) function_2d_name(j), center(1:dim_num), &
        radius1, radius2, result1, result3 - result2

    end do

  end do

  return
end
subroutine test085

!*****************************************************************************80
!
!! TEST085 tests CIRCLE_ANNULUS_AREA_2D.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: dim_num = 2
  integer ( kind = 4 ), parameter :: test_num = 3

  real    ( kind = 8 ) area
  real    ( kind = 8 ) center(dim_num)
  real    ( kind = 8 ), dimension ( dim_num, test_num ) :: center_test = reshape ( (/ &
    0.0D+00, 0.0D+00, &
    1.0D+00, 0.0D+00, &
    3.0D+00, 4.0D+00 /), (/ dim_num, test_num /) )
  real    ( kind = 8 ) circle_annulus_area_2d
  integer ( kind = 4 ) i
  real    ( kind = 8 ) radius1
  real    ( kind = 8 ), dimension ( test_num ) :: radius1_test = (/ &
    0.0D+00, 1.0D+00, 1.0D+00 /)
  real    ( kind = 8 ) radius2
  real    ( kind = 8 ), dimension ( test_num ) :: radius2_test = (/ &
    1.0D+00, 2.0D+00, 3.0D+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST085'
  write ( *, '(a)' ) '  CIRCLE_ANNULUS_AREA_2D computes the area of a '
  write ( *, '(a)' ) '    circular annulus.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '      CENTER       Radius1   Radius2   Area'
  write ( *, '(a)' ) ' '

  do i = 1, test_num

    center(1:dim_num) = center_test(1:dim_num,i)

    radius1 = radius1_test(i)
    radius2 = radius2_test(i)

    area = circle_annulus_area_2d ( radius1, radius2 )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,5f11.5)' ) center(1:dim_num), radius1, radius2, area

  end do

  return
end
subroutine test09

!*****************************************************************************80
!
!! TEST09 tests CIRCLE_ANNULUS_SECTOR, CIRCLE_RT_SET, CIRCLE_RT_SUM.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2
  integer   ( kind = 4 ), parameter :: test_num = 4

  real      ( kind = 8 ) as1
  real      ( kind = 8 ) as2
  real      ( kind = 8 ) as3
  real      ( kind = 8 ) as4
  real      ( kind = 8 ) center(dim_num)
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) nc
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) nr
  integer   ( kind = 4 ) nr2
  integer   ( kind = 4 ) nt
  real      ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real      ( kind = 8 ) ra(5)
  real      ( kind = 8 ) radius
  real      ( kind = 8 ) radius1a
  real      ( kind = 8 ) radius2a
  real      ( kind = 8 ) radius1b
  real      ( kind = 8 ) radius2b
  real      ( kind = 8 ) radius1c
  real      ( kind = 8 ) radius2c
  real      ( kind = 8 ) radius1d
  real      ( kind = 8 ) radius2d
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ) rw(5)
  real      ( kind = 8 ) ta(20)
  real      ( kind = 8 ) theta1a
  real      ( kind = 8 ) theta2a
  real      ( kind = 8 ) theta1b
  real      ( kind = 8 ) theta2b
  real      ( kind = 8 ) theta1c
  real      ( kind = 8 ) theta2c
  real      ( kind = 8 ) theta1d
  real      ( kind = 8 ) theta2d
  real      ( kind = 8 ) tw(20)
  real      ( kind = 8 ) zw

  nr = 5

  rule = 9
  call circle_rt_size ( rule, nr2, nt, nc )
  call circle_rt_set ( rule, nr2, nt, nc, ra, rw, ta, tw, zw )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  CIRCLE_ANNULUS_SECTOR estimates an integral in a '
  write ( *, '(a)' ) '    circular annulus sector.'
  write ( *, '(a)' ) '  CIRCLE_RT_SET sets an integration rule in a circle.'
  write ( *, '(a)' ) '  CIRCLE_RT_SUM uses an integration rule in a circle.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  To test CIRCLE_ANNULUS_SECTOR, we estimate an integral'
  write ( *, '(a)' ) '  over 4 annular sectors that make up the unit circle, '
  write ( *, '(a)' ) '  and add to get RESULT1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We will also estimate the integral over the unit circle'
  write ( *, '(a)' ) '  using CIRCLE_RT_SET and CIRCLE_RT_SUM to get RESULT2.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We will then compare RESULT1 and RESULT2.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  CIRCLE_ANNULUS_SECTOR computations will use NR = ',nr
  write ( *, '(a,i8)' ) '  CIRCLE_RT_SET/CIRCLE_RT_SUM will use rule ', rule
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  "RESULT1" is the sum of Annulus Sector calculations.'
  write ( *, '(a)' ) '  "RESULT2" is for CIRCLE_RT_SET/CIRCLE_RT_SUM.'
  write ( *, '(a)' ) ' '

  center(1:dim_num) = (/ 0.0D+00, 0.0D+00 /)
  radius = 1.0D+00

  radius1a = 0.0D+00
  radius2a = 0.25D+00
  theta1a = 0.0D+00
  theta2a = 0.5D+00 * pi

  radius1b = 0.0D+00
  radius2b = 0.25D+00
  theta1b = 0.5D+00 * pi
  theta2b = 2.0D+00 * pi

  radius1c = 0.25D+00
  radius2c = 1.0D+00
  theta1c = 0.0D+00
  theta2c = 0.25D+00 * pi

  radius1d = 0.25D+00
  radius2d = 1.0D+00
  theta1d = 0.25D+00 * pi
  theta2d = 2.0D+00 * pi

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '       F  Result1  Result2'
  write ( *, '(a)' ) ' '

  num = function_2d_num ( )

  do j = 1, num

    call function_2d_set ( 'SET', j )

    call circle_annulus_sector ( function_2d, center, radius1a, radius2a, theta1a, &
      theta2a, nr, as1 )

    call circle_annulus_sector ( function_2d, center, radius1b, radius2b, theta1b, &
      theta2b, nr, as2 )

    call circle_annulus_sector ( function_2d, center, radius1c, radius2c, theta1c, &
      theta2c, nr, as3 )

    call circle_annulus_sector ( function_2d, center, radius1d, radius2d, theta1d, &
      theta2d, nr, as4 )

    result1 = as1 + as2 + as3 + as4

    call circle_rt_sum ( function_2d, center, radius, nr2, ra, rw, nt, ta, tw, zw, &
      result2 )

    write ( *, '(2x,a7,2g14.6)' ) function_2d_name(j), result1, result2

  end do

  return
end
subroutine test10

!*****************************************************************************80
!
!! TEST10 tests CIRCLE_CUM.
!
!  Modified:
!
!    07 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2

  real      ( kind = 8 ), dimension ( dim_num ) :: center = (/ 0.0D+00, 0.0D+00 /)
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) :: r = 3.0D+00
  real      ( kind = 8 ) result(4)
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  CIRCLE_CUM approximates an integral over a circle.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  We use radius R = ', r
  write ( *, '(a)' ) '  and center:'
  write ( *, '(a,2g14.6)' ) '  CENTER = ', center(1:dim_num)
  write ( *, '(a)' ) ' '

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    Order:      2             4              8     ' // &
    '       16'
  write ( *, '(a)' ) '  F(X)'
  write ( *, '(a)' ) ' '
 
  num = function_2d_num ( )

  do i = 1, num

    call function_2d_set ( 'SET', i )

    do j = 1, 4

      order = 2**j

      call circle_cum ( function_2d, center, r, order, result(j) )

    end do

    write ( *, '(2x,a7,5f14.8)' ) function_2d_name(i), result(1:4)
 
  end do

  return
end
subroutine test11

!*****************************************************************************80
!
!! TEST11 tests LENS_HALF_AREA_2D, CIRCLE_SECTOR_AREA_2D, CIRCLE_TRIANGLE_AREA_2D.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real    ( kind = 8 ) area1
  real    ( kind = 8 ) area2
  real    ( kind = 8 ) area3
  real    ( kind = 8 ) lens_half_area_2d
  real    ( kind = 8 ) circle_sector_area_2d
  real    ( kind = 8 ) circle_triangle_area_2d
  integer ( kind = 4 ) i
  real    ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real    ( kind = 8 ) r
  real    ( kind = 8 ) theta1
  real    ( kind = 8 ) theta2

  r = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  LENS_HALF_AREA_2D computes the area of a'
  write ( *, '(a)' ) '    circular half lens, defined by joining the endpoints'
  write ( *, '(a)' ) '    of a circular arc.'
  write ( *, '(a)' ) '  CIRCLE_SECTOR_AREA_2D computes the area of a'
  write ( *, '(a)' ) '    circular sector, defined by joining the endpoints'
  write ( *, '(a)' ) '    of a circular arc to the center.'
  write ( *, '(a)' ) '  CIRCLE_TRIANGLE_AREA_2D computes the signed area of a'
  write ( *, '(a)' ) '    triangle, defined by joining the endpoints'
  write ( *, '(a)' ) '    of a circular arc and the center.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '      R       Theta1 Theta2        ' // &
    'Sector       Triangle     Half Lens'
  write ( *, '(a)' ) ' '

  do i = 0, 12

    theta1 = 0.0D+00
    theta2 = real ( i, kind = 8 ) * 2.0D+00 * pi / 12.0D+00

    area1 = circle_sector_area_2d ( r, theta1, theta2 )

    area2 = circle_triangle_area_2d ( r, theta1, theta2 )

    area3 = lens_half_area_2d ( r, theta1, theta2 )

    write ( *, '(2x,f9.3,f9.3,4f14.8)' ) r, theta1, theta2, area1, area2, area3

  end do

  return
end
subroutine test12

!*****************************************************************************80
!
!! TEST12 tests LENS_HALF_AREA_2D, LENS_HALF_H_AREA_2D, LENS_HALF_W_AREA_2D.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real    ( kind = 8 ) area1
  real    ( kind = 8 ) area2
  real    ( kind = 8 ) area3
  real    ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real    ( kind = 8 ) h
  integer ( kind = 4 ) i
  real    ( kind = 8 ) lens_half_area_2d
  real    ( kind = 8 ) lens_half_h_area_2d
  real    ( kind = 8 ) lens_half_w_area_2d
  real    ( kind = 8 ) r
  real    ( kind = 8 ) theta1
  real    ( kind = 8 ) theta2
  real    ( kind = 8 ) w

  r = 50.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  For the area of a circular half lens,'
  write ( *, '(a)' ) '  LENS_HALF_AREA_2D uses two angles;'
  write ( *, '(a)' ) '  LENS_HALF_H_AREA_2D works from the height;'
  write ( *, '(a)' ) '  LENS_HALF_W_AREA_2D works from the width.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  The circle has radius R = ', r
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  THETA1 THETA2  H     W  Area(THETA) Area(H)  Area(W)'
  write ( *, '(a)' ) ' '

  do i = 0, 12

    theta1 = 0.0D+00
    theta2 = real ( i, kind = 8 ) * 2.0D+00 * pi / 12.0D+00
    w = 2.0D+00 * r * sin ( 0.5D+00 * ( theta2 - theta1 ) ) 
    h = r * ( 1.0D+00 - cos ( 0.5D+00 * ( theta2 - theta1 ) ) )

    area1 = lens_half_area_2d ( r, theta1, theta2 )

    area2 = lens_half_h_area_2d ( r, h )

    area3 = lens_half_w_area_2d ( r, w )

    write ( *, '(2x,4f6.2,3f10.4)' ) theta1, theta2, h, w, area1, area2, area3

  end do

  return
end
subroutine test13

!*****************************************************************************80
!
!! TEST13 tests CIRCLE_SECTOR, CIRCLE_SECTOR_AREA_2D.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2
  integer   ( kind = 4 ), parameter :: nrlo = 1
  integer   ( kind = 4 ), parameter :: nrhi = 5
  integer   ( kind = 4 ), parameter :: test_num = 4

  real      ( kind = 8 ) area
  real      ( kind = 8 ) center(dim_num)
  real      ( kind = 8 ), dimension ( dim_num, test_num ) :: center_test = reshape ( (/ &
    0.0D+00, 0.0D+00, &
    0.0D+00, 0.0D+00, &
    0.0D+00, 0.0D+00, &
    0.0D+00, 0.0D+00 /), (/ dim_num, test_num /) )
  real      ( kind = 8 ) circle_sector_area_2d
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) nr
  real      ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real      ( kind = 8 ) radius
  real      ( kind = 8 ), dimension ( test_num ) :: radius_test = (/ &
    1.0D+00, 2.0D+00, 4.0D+00, 8.0D+00 /)
  real      ( kind = 8 ) result(nrlo:nrhi)
  real      ( kind = 8 ) theta1
  real      ( kind = 8 ), dimension ( test_num ) :: theta1_test = (/ &
    0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 /)
  real      ( kind = 8 ) theta2
  real      ( kind = 8 ), dimension ( test_num ) :: theta2_test = (/ &
    2.0D+00, 1.0D+00, 0.5D+00, 0.25D+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST13'
  write ( *, '(a)' ) '  CIRCLE_SECTOR_AREA_2D computes the area '
  write ( *, '(a)' ) '    of a circular sector.'
  write ( *, '(a)' ) '  CIRCLE_SECTOR estimates an integral '
  write ( *, '(a)' ) '    in a circular sector.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The user can specify NR, the number of radial values'
  write ( *, '(a)' ) '  used to approximated the integral.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, computations will use values of NR'
  write ( *, '(a,i8)' ) '  from ', nrlo
  write ( *, '(a,i8)' ) '  to   ', nrhi
  write ( *, '(a)' ) ' '

  do i = 1, test_num

    center(1:dim_num) = center_test(1:dim_num,i)

    radius = radius_test(i)
    theta1 = theta1_test(i) * pi
    theta2 = theta2_test(i) * pi

    area = circle_sector_area_2d ( radius, theta1, theta2 )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '       CENTER      RADIUS  THETA1  THETA2  Area'
    write ( *, '(a)' ) ' '
    write ( *, '(6f8.4)' ) center(1:dim_num), radius, theta1, theta2, area
    write ( *, '(a)' ) ' '
    write ( *, '(a7,14(6x,i2,6x))' ) '   F   ', ( nr, nr = nrlo, nrhi )
    write ( *, '(a)' ) ' '

    num = function_2d_num ( )

    do j = 1, num

      call function_2d_set ( 'SET', j )

      do nr = nrlo, nrhi
        call circle_sector ( function_2d, center, radius, theta1, theta2, nr, &
          result(nr) )
      end do

      write ( *, '(2x,a7,5g14.6)' ) function_2d_name(j), result(nrlo:nrhi)

    end do

  end do

  return
end
subroutine test14

!*****************************************************************************80
!
!! TEST14 tests CIRCLE_SECTOR, CIRCLE_RT_SET, CIRCLE_RT_SUM.
!
!  Modified:
!
!    07 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2
  integer   ( kind = 4 ), parameter :: test_num = 4

  real      ( kind = 8 ) area1
  real      ( kind = 8 ) area2
  real      ( kind = 8 ) area3
  real      ( kind = 8 ) center(dim_num)
  real      ( kind = 8 ), dimension ( dim_num, test_num ) :: center_test = reshape ( (/ &
    0.0D+00, 0.0D+00, &
    0.0D+00, 0.0D+00, &
    0.0D+00, 0.0D+00, &
    0.0D+00, 0.0D+00 /), (/ dim_num, test_num /) )
  real      ( kind = 8 ) circle_area_2d
  real      ( kind = 8 ) circle_sector_area_2d
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) nc
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) nr
  integer   ( kind = 4 ) nr2
  integer   ( kind = 4 ) nt
  real      ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real      ( kind = 8 ) ra(5)
  real      ( kind = 8 ) radius
  real      ( kind = 8 ), dimension ( test_num ) :: radius_test = (/ &
    1.0D+00, 2.0D+00, 4.0D+00, 8.0D+00 /)
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) resulta
  real      ( kind = 8 ) resultb
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ) rw(5)
  real      ( kind = 8 ) ta(20)
  real      ( kind = 8 ) theta1
  real      ( kind = 8 ), dimension ( test_num ) :: theta1_test = (/ &
    0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 /)
  real      ( kind = 8 ) theta2
  real      ( kind = 8 ), dimension ( test_num ) :: theta2_test = (/ &
    2.0D+00, 1.0D+00, 0.5D+00, 0.25D+00 /)
  real      ( kind = 8 ) theta3
  real      ( kind = 8 ) tw(20)
  real      ( kind = 8 ) zw

  nr = 5

  rule = 9
  call circle_rt_size ( rule, nr2, nt, nc )
  call circle_rt_set ( rule, nr2, nt, nc, ra, rw, ta, tw, zw )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST14'
  write ( *, '(a)' ) '  CIRCLE_SECTOR estimates integrals in a circular sector.'
  write ( *, '(a)' ) '  CIRCLE_RT_SET sets an integration rule in a circle.'
  write ( *, '(a)' ) '  CIRCLE_RT_SUM uses an integration rule in a circle.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  To test CIRCLE_SECTOR, we estimate an integral over'
  write ( *, '(a)' ) '  a sector, and over its complement and add the results'
  write ( *, '(a)' ) '  to get RESULT1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We also estimate the integral over the whole circle'
  write ( *, '(a)' ) '  using CIRCLE_RT_SET and CIRCLE_RT_SUM to get RESULT2.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We will then compare RESULT1 and RESULT2.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  CIRCLE_SECTOR computations will use NR = ',nr
  write ( *, '(a,i8)' ) '  CIRCLE_RT_SET/CIRCLE_RT_SUM will use rule ', rule
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  "Sector1" and "Sector2" are the CIRCLE_SECTOR'
  write ( *, '(a)' ) '  computations'
  write ( *, '(a)' ) '  for the sector and its complement.'
  write ( *, '(a)' ) '  "Sum" is the sum of Sector1 and Sector2.'
  write ( *, '(a)' ) '  "Circle" is the computation for '
  write ( *, '(a)' ) '  CIRCLE_RT_SET + CIRCLE_RT_SUM.'
  write ( *, '(a)' ) ' '

  do i = 1, test_num

    center(1:dim_num) = center_test(1:dim_num,i)

    radius = radius_test(i)

    theta1 = theta1_test(i) * pi
    theta2 = theta2_test(i) * pi
    theta3 = theta2 + 2.0D+00 * pi - ( theta2 - theta1 )

    area1 = circle_sector_area_2d ( radius, theta1, theta2 )
    area2 = circle_sector_area_2d ( radius, theta2, theta3 )
    area3 = circle_area_2d ( radius )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '      CENTER       RADIUS   THETA1   THETA2   Area1   Area2  Circle'
    write ( *, '(a)' ) ' '
    write ( *, '(8f9.4)' ) center(1:dim_num), radius, theta1, theta2, area1, area2, area3
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '       F   Sector1       Sector2         Sum         Circle'
    write ( *, '(a)' ) ' '

    num = function_2d_num ( )

    do j = 1, num

      call function_2d_set ( 'SET', j )

      call circle_sector ( function_2d, center, radius, theta1, theta2, nr, &
        resulta )

      call circle_sector ( function_2d, center, radius, theta2, theta3, nr, &
        resultb )

      result1 = resulta + resultb

      call circle_rt_sum ( function_2d, center, radius, nr2, ra, rw, nt, ta, &
        tw, zw, result2 )

      write ( *, '(2x,a7,4g14.6)' ) function_2d_name(j), resulta, resultb, &
        result1, result2

    end do

  end do

  return
end
subroutine test15

!*****************************************************************************80
!
!! TEST15 tests CIRCLE_RT_SET and CIRCLE_RT_SUM.
!
!  Modified:
!
!    19 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2
  integer   ( kind = 4 ), parameter :: rule_max = 9

  real      ( kind = 8 ), dimension ( dim_num ) :: center = (/ 1.0D+00, 1.0D+00 /)
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) nc
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) nr
  integer   ( kind = 4 ) nt
  real      ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real      ( kind = 8 ), parameter :: r = 1.0D+00
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ra
  real      ( kind = 8 ) result(rule_max)
  real      ( kind = 8 ), allocatable, dimension ( : ) :: rw
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ta
  real      ( kind = 8 ), allocatable, dimension ( : ) :: tw
  real      ( kind = 8 ) zw

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST15'
  write ( *, '(a)' ) '  For R, Theta product rules on the unit circle,'
  write ( *, '(a)' ) '  CIRCLE_RT_SET sets a rule.'
  write ( *, '(a)' ) '  CIRCLE_RT_SUM uses the rule in an arbitrary circle.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  We use a radius ', r
  write ( *, '(a)' ) '  and center:'
  write ( *, '(a,2g14.6)' ) '  CENTER = ', center(1:dim_num)
  write ( *, '(a)' ) ' '

  do ilo = 1, rule_max, 5

    ihi = min ( ilo +  4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a7,7x,5(i7,7x))' ) 'Rule:  ', ( rule, rule = ilo, ihi )
    write ( *, '(2x,a)' ) 'Function '

    num = function_2d_num ( )

    do i = 1, num

      call function_2d_set ( 'SET', i )

      do rule = ilo, ihi

        call circle_rt_size ( rule, nr, nt, nc )

        allocate ( ra(1:nr) )
        allocate ( rw(1:nr) )
        allocate ( ta(1:nt) )
        allocate ( tw(1:nt) )

        call circle_rt_set ( rule, nr, nt, nc, ra, rw, ta, tw, zw )

        call circle_rt_sum ( function_2d, center, r, nr, ra, rw, nt, ta, tw, zw, &
          result(rule) )

        deallocate ( ra )
        deallocate ( rw )
        deallocate ( ta )
        deallocate ( tw )

      end do

      write ( *, '(2x,a7,5f14.8)' ) function_2d_name(i), result(ilo:ihi)

    end do

  end do

  return
end
subroutine test16

!*****************************************************************************80
!
!! TEST16 tests CIRCLE_XY_SET and CIRCLE_XY_SUM.
!
!  Modified:
!
!    16 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2
  integer   ( kind = 4 ), parameter :: rule_max = 13

  real      ( kind = 8 ), dimension ( dim_num ) :: center = (/ 1.0D+00, 1.0D+00 /)
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) r
  real      ( kind = 8 ) result(rule_max)
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab

  r = 1.0D+00
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST16'
  write ( *, '(a)' ) '  CIRCLE_XY_SET sets a quadrature rule '
  write ( *, '(a)' ) '    for the unit circle.'
  write ( *, '(a)' ) '  CIRCLE_XY_SUM evaluates the quadrature rule'
  write ( *, '(a)' ) '    in an arbitrary circle.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  We use a radius ', r
  write ( *, '(a)' ) '  and center:'
  write ( *, '(a,2g14.6)' ) '  CENTER = ', center(1:dim_num)
  write ( *, '(a)' ) ' '

  do ilo = 1, rule_max, 5

    ihi = min ( ilo +  4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a7,7x,5(i7,7x))' ) 'Rule:  ', ( rule, rule = ilo, ihi )
    write ( *, '(2x,a)' ) 'Function '

    num = function_2d_num ( )

    do i = 1, num

      call function_2d_set ( 'SET', i )

      do rule = ilo, ihi

        call circle_xy_size ( rule, order )

        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( weight(1:order) )

        call circle_xy_set ( rule, order, xtab, ytab, weight )

        call circle_xy_sum ( function_2d, center, r, order, xtab, ytab, weight, &
          result(rule) )

        deallocate ( weight )
        deallocate ( xtab )
        deallocate ( ytab )

      end do

      write ( *, '(2x,a7,5f14.8)' ) function_2d_name(i), result(ilo:ihi)

    end do

  end do

  return
end
subroutine test17

!*****************************************************************************80
!
!! TEST17 tests CONE_UNIT_3D, CONE_VOLUME_3D.
!
!  Modified:
!
!    05 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ) cone_volume_3d
  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  real      ( kind = 8 ) h
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) num
  real      ( kind = 8 ), parameter :: r = 1.0D+00
  real      ( kind = 8 ) result

  h = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST17'
  write ( *, '(a)' ) '  CONE_UNIT_3D approximates integrals in a unit cone.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Volume = ', cone_volume_3d ( r, h )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    F(X)    CONE_3D'
  write ( *, '(a)' ) ' '

  num = function_3d_num ( )

  do i = 1, num

    call function_3d_set ( 'SET', i )

    call cone_unit_3d ( function_3d, result )

    write ( *, '(2x,a7,f14.8)' ) function_3d_name(i), result

  end do

  return
end
subroutine test18

!*****************************************************************************80
!
!! TEST18 tests CUBE_SHELL_ND.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: n_max = 4
  integer   ( kind = 4 ), parameter :: test_num = 2

  real      ( kind = 8 ) cube_shell_volume_nd
  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) function_nd_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) r1
  real      ( kind = 8 ), dimension ( test_num ) :: r1_test = (/ &
    0.0D+00, 1.0D+00 /)
  real      ( kind = 8 ) :: r2
  real      ( kind = 8 ), dimension ( test_num ) :: r2_test = (/ &
    1.0D+00, 2.0D+00 /)
  real      ( kind = 8 ) result
  integer   ( kind = 4 ) test

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST18'
  write ( *, '(a)' ) '  CUBE_SHELL_ND approximates integrals in a'
  write ( *, '(a)' ) '    cubical shell in ND.'

  do test = 1, test_num

    r1 = r1_test(test)
    r2 = r2_test(test)

    write ( *, '(a)' ) ' '
    write ( *, '(a,g14.6)' ) '  Inner radius = ', r1
    write ( *, '(a,g14.6)' ) '  Outer radius = ', r2
    write ( *, '(a)' ) ' '
 
    do n = 2, n_max
 
      write ( *, '(a)' ) ' '
      write ( *, '(a,i8)' ) '  Spatial dimension N = ', n
      write ( *, '(a,g14.6)' ) '  Volume = ', cube_shell_volume_nd ( n, r1, r2 )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '    F(X)      CUBE_SHELL_ND'
      write ( *, '(a)' ) ' '
 
      num = function_nd_num ( )

      do i = 1, num

        call function_nd_set ( 'SET', i )

        call cube_shell_nd ( function_nd, n, r1, r2, result )

        write ( *, '(2x,a7,f14.8)' ) function_nd_name(i), result

      end do

    end do
 
  end do
 
  return
end
subroutine test19

!*****************************************************************************80
!
!! TEST19 tests CUBE_UNIT_3D, QMULT_3D, RECTANGLE_3D.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ) a1
  real      ( kind = 8 ) a(3)
  real      ( kind = 8 ) b1
  real      ( kind = 8 ) b(3)
  real      ( kind = 8 ), external :: fl18
  real      ( kind = 8 ), external :: fl28
  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: fu18
  real      ( kind = 8 ), external :: fu28
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ), parameter :: n = 3
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) qmult_3d
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) result3

  a1 = -1.0D+00
  b1 = +1.0D+00

  a(1) = -1.0D+00
  a(2) = -1.0D+00
  a(3) = -1.0D+00
  b(1) = 1.0D+00
  b(2) = 1.0D+00
  b(3) = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST19'
  write ( *, '(a)' ) '  CUBE_UNIT_3D approximates integrals '
  write ( *, '(a)' ) '    in the unit cube in 3D.'
  write ( *, '(a)' ) '  QMULT_3D approximates triple integrals.'
  write ( *, '(a)' ) '  RECTANGLE_3D approximates integrals '
  write ( *, '(a)' ) '    in a rectangular block.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    F(X)    CUBE_UNIT_3D  QMULT_3D      RECTANGLE_3D'
  write ( *, '(a)' ) ' '

  num = function_3d_num ( )

  do i = 1, num

    call function_3d_set ( 'SET', i )

    call cube_unit_3d ( function_3d, result1 )
    result2 = qmult_3d ( function_3d, a1, b1, fu18, fl18, fu28, fl28 )
    call rectangle_3d ( function_3d, a, b, result3 )

    write ( *, '(2x,a7,3f14.8)' ) function_3d_name(i), result1, result2, result3

  end do

  return
end
subroutine test20

!*****************************************************************************80
!
!! TEST20 tests CUBE_UNIT_ND.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: max_k = 10
  integer   ( kind = 4 ), parameter :: max_test = 2

  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) function_nd_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) i_test
  integer   ( kind = 4 ) k
  integer   ( kind = 4 ) khi
  integer   ( kind = 4 ) klo
  integer   ( kind = 4 ), parameter, dimension ( max_test ) :: k_test = (/ 10, 5 /)
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ), parameter, dimension ( max_test ) :: n_test = (/ 2, 3 /)
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) qa(max_k)
  real      ( kind = 8 ) qb(max_k)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST20'
  write ( *, '(a)' ) '  CUBE_UNIT_ND approximates integrals inside '
  write ( *, '(a)' ) '    the unit cube in ND.'

  do i_test = 1, max_test
 
    n = n_test(i_test)
    k = k_test(i_test)
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Spatial dimension N = ', n
    write ( *, '(a,i8)' ) '  Value of K = ', k
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '    F(X)    CUBE_UNIT_ND'
    write ( *, '(a)' ) ' '

    num = function_nd_num ( )

    do i = 1, num

      call function_nd_set ( 'SET', i )

      call cube_unit_nd ( function_nd, qa, qb, n, k )

      do klo = 1, k, 5
        khi = min ( klo + 4, k )
        if ( klo == 1 ) then
          write ( *, '(2x,a7,5f14.8)' ) function_nd_name(i), qa(klo:khi)
        else
          write ( *, '(2x,7x,5f14.8)' )           qa(klo:khi)
        end if
      end do

      do klo = 1, k, 5
        khi = min ( klo + 4, k )
        write ( *, '(2x,7x,5f14.8)' )           qb(klo:khi)
      end do
                                           
    end do

  end do

  return
end
subroutine test205

!*****************************************************************************80
!
!! TEST205 tests ELLIPSE_AREA_2D, ELLIPSE_CIRCUMFERENCE_2D, ELLIPSE_ECCENTRICITY_2D.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real    ( kind = 8 ) area
  real    ( kind = 8 ) r8_uniform_01
  real    ( kind = 8 ) e
  real    ( kind = 8 ) ellipse_area_2d
  real    ( kind = 8 ) ellipse_circumference_2d
  real    ( kind = 8 ) ellipse_eccentricity_2d
  integer ( kind = 4 ) i
  real    ( kind = 8 ) p
  real    ( kind = 8 ) r1
  real    ( kind = 8 ) r2
  integer ( kind = 4 ) :: seed = 123456789

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST205'
  write ( *, '(a)' ) '  ELLIPSE_AREA_2D returns the area of an ellipse.'
  write ( *, '(a)' ) '  ELLIPSE_ECCENTRICITY_2D returns the '
  write ( *, '(a)' ) '    eccentricity of an ellipse.'
  write ( *, '(a)' ) '  ELLIPSE_CIRCUMFERENCE_2D returns the '
  write ( *, '(a)' ) '    circumference of an ellipse.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '        R1        R2        E         Circum    Area'
  write ( *, '(a)' ) ' '

  do i = 1, 5

    if ( i == 1 ) then
      r1 = 25.0D+00
      r2 = 20.0D+00
    else
      r1 = r8_uniform_01 ( seed )
      r2 = r8_uniform_01 ( seed )
    end if

    e = ellipse_eccentricity_2d ( r1, r2 )
    p = ellipse_circumference_2d ( r1, r2 )    
    area = ellipse_area_2d ( r1, r2 )

    write ( *, '(2x,5f10.4)' ) r1, r2, e, p, area

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  (For the first example, '
  write ( *, '(a)' ) '  the eccentricity should be 0.6,'
  write ( *, '(a)' ) '  the circumference should be about 141.8).'

  return
end
subroutine test21

!*****************************************************************************80
!
!! TEST21 tests HEXAGON_UNIT_SET and HEXAGON_SUM.
!
!  Modified:
!
!    13 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2
  integer   ( kind = 4 ), parameter :: rule_max = 4

  real      ( kind = 8 ), dimension ( dim_num ) :: center = (/ &
    0.0D+00, 0.0D+00 /)
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) rad
  real      ( kind = 8 ) result(rule_max)
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab

  rad = 2.0D+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST21'
  write ( *, '(a)' ) '  HEXAGON_UNIT_SET sets a quadrature rule for the '
  write ( *, '(a)' ) '    unit hexagon.'
  write ( *, '(a)' ) '  HEXAGON_SUM evaluates the quadrature rule'
  write ( *, '(a)' ) '    in an arbitrary hexagon.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  We use a radius ', rad
  write ( *, '(a)' ) '  and center:'
  write ( *, '(a,2g14.6)' ) '  CENTER = ', center(1:dim_num)
  write ( *, '(a)' ) ' '

  do ilo = 1, rule_max, 5

    ihi = min ( ilo + 4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a,5i6)' ) 'Rule:    ', ( rule, rule = ilo, ihi )
    write ( *, '(2x,a)' ) 'Function '
    write ( *, '(a)' ) ' '

    num = function_2d_num ( )

    do i = 1, num

      call function_2d_set ( 'SET', i )

      do rule = ilo, ihi

        call hexagon_unit_size ( rule, order )

        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( weight(1:order) )

        call hexagon_unit_set ( rule, order, xtab, ytab, weight )

        call hexagon_sum ( function_2d, center, rad, order, xtab, ytab, weight, &
          result(rule) )

        deallocate ( xtab )
        deallocate ( ytab )
        deallocate ( weight )

      end do

      write ( *, '(2x,a7,5f14.8)' ) function_2d_name(i), result(ilo:ihi)

    end do

  end do

  return
end
subroutine test215

!*****************************************************************************80
!
!! TEST215 tests LENS_HALF_2D.
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: dim_num = 2

  real    ( kind = 8 ) area
  real    ( kind = 8 ), dimension ( dim_num ) :: center = (/ &
    0.0D+00, 0.0D+00 /)
  real    ( kind = 8 ), external :: f_1_2d
  real    ( kind = 8 ), external :: f_x_2d
  real    ( kind = 8 ), external :: f_r_2d
  integer ( kind = 4 ) i
  integer ( kind = 4 ) :: i_max = 8
  real    ( kind = 8 ) lens_half_2d
  real    ( kind = 8 ) lens_half_area_2d
  integer ( kind = 4 ) order
  real    ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real    ( kind = 8 ) r
  real    ( kind = 8 ) theta1
  real    ( kind = 8 ) theta2
  real    ( kind = 8 ) value

  r = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST215'
  write ( *, '(a)' ) '  LENS_HALF_2D approximates an integral within a'
  write ( *, '(a)' ) '    circular half lens, defined by joining the endpoints'
  write ( *, '(a)' ) '    of a circular arc.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integrate F(X,Y) = 1'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '      R            Theta1      Theta2        ' // &
    'Area        Order Integral'
  write ( *, '(a)' ) ' '

  do i = 0, i_max

    theta1 = 0.0D+00
    theta2 = real ( i, kind = 8 ) * 2.0D+00 * pi &
      / real ( i_max, kind = 8 )

    area = lens_half_area_2d ( r, theta1, theta2 )

    write ( *, '(a)' ) ' '

    do order = 2, 16, 2
      value = lens_half_2d ( f_1_2d, center, r, theta1, theta2, order )
      write ( *, '(4g14.6,i8,g14.6)' ) r, theta1, theta2, area, order, value
    end do

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integrate F(X,Y) = X'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '      R            Theta1      Theta2        ' // &
    'Area        Order Integral'
  write ( *, '(a)' ) ' '

  do i = 0, i_max

    theta1 = 0.0D+00
    theta2 = real ( i, kind = 8 ) * 2.0D+00 * pi &
      / real ( i_max, kind = 8 )

    area = lens_half_area_2d ( r, theta1, theta2 )

    write ( *, '(a)' ) ' '

    do order = 2, 16, 2
      value = lens_half_2d ( f_x_2d, center, r, theta1, theta2, order )
      write ( *, '(4g14.6,i8,g14.6)' ) r, theta1, theta2, area, order, value
    end do

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integrate F(X,Y) = R'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '      R            Theta1      Theta2        ' // &
    'Area        Order Integral'
  write ( *, '(a)' ) ' '

  do i = 0, i_max

    theta1 = 0.0D+00
    theta2 = real ( i, kind = 8 ) * 2.0D+00 * pi &
      / real ( i_max, kind = 8 )

    area = lens_half_area_2d ( r, theta1, theta2 )

    write ( *, '(a)' ) ' '

    do order = 2, 16, 2
      value = lens_half_2d ( f_r_2d, center, r, theta1, theta2, order )
      write ( *, '(4g14.6,i8,g14.6)' ) r, theta1, theta2, area, order, value
    end do

  end do

  return
end
subroutine test22

!*****************************************************************************80
!
!! TEST22 tests OCTAHEDRON_UNIT_ND.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: n_max = 3

  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) function_nd_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result(3)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST22'
  write ( *, '(a)' ) '  OCTAHEDRON_UNIT_ND approximates integrals in a unit' 
  write ( *, '(a)' ) '    octahedron in N dimensions.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    F(X)    N = 1    N = 2   N = 3 '
  write ( *, '(a)' ) ' '
 
  num = function_nd_num ( )

  do i = 1, num

    call function_nd_set ( 'SET', i )

    do n = 1, n_max
      call octahedron_unit_nd ( function_nd, n, result(n) )
    end do

    write ( *, '(2x,a7,3f14.8)' ) function_nd_name(i), result(1:n_max)

  end do

  return
end
subroutine test23

!*****************************************************************************80
!
!! TEST23 tests PARALLELIPIPED_VOLUME_ND.
!
!  Modified:
!
!    06 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) i
  integer ( kind = 4 ) n
  real    ( kind = 8 ) parallelipiped_volume_nd
  real    ( kind = 8 ), allocatable, dimension ( :, : ) :: v
  real    ( kind = 8 ) volume

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST23'
  write ( *, '(a)' ) '  PARALLELIPIPED_VOLUME_ND computes the volume of a'
  write ( *, '(a)' ) '    parallelipiped in N dimensions.'
  write ( *, '(a)' ) ' '

  do n = 2, 4

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Spatial dimension N = ', n

    allocate ( v(1:n,1:n+1) )
!
!  Set the values of the parallelipiped.
!
    call setsim ( n, v )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Parallelipiped vertices:'
    write ( *, '(a)' ) ' '

    do i = 1, n
      write ( *, '(2x,5f4.0)' ) v(i,1:n+1)
    end do

    volume = parallelipiped_volume_nd ( n, v )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a,g14.6)' ) '  Volume is ', volume

    deallocate ( v )

  end do

  return
end
subroutine test24

!*****************************************************************************80
!
!! TEST24 tests POLYGON_**_2D;
!
!  Modified:
!
!    10 December 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: npts = 4

  character ( len = 4 ) name
  real      ( kind = 8 ) result
  real      ( kind = 8 ), dimension ( npts) :: x = (/ &
    0.0D+00, 1.0D+00, 1.0D+00, 0.0D+00 /)
  real      ( kind = 8 ), dimension ( npts) :: y = (/ &
    0.0D+00, 0.0D+00, 1.0D+00, 1.0D+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST24'
  write ( *, '(a)' ) '  For a polygon in 2D:' 
  write ( *, '(a)' ) '  POLYGON_1_2D integrates 1'
  write ( *, '(a)' ) '  POLYGON_X_2D integrates X'
  write ( *, '(a)' ) '  POLYGON_Y_2D integrates Y'
  write ( *, '(a)' ) '  POLYGON_XX_2D integrates X**2'
  write ( *, '(a)' ) '  POLYGON_XY_2D integrates X*Y'
  write ( *, '(a)' ) '  POLYGON_YY_2D integrates Y**2'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  F(X,Y)    Integral'
  write ( *, '(a)' ) ' '
 
  name = '1'
  call polygon_1_2d ( npts, x, y, result )
  write ( *, '(2x,a4,4x,g14.6)' ) name, result
 
  name = 'X'
  call polygon_x_2d ( npts, x, y, result )
  write ( *, '(2x,a4,4x,g14.6)' ) name, result
 
  name = 'Y'
  call polygon_y_2d ( npts, x, y, result )
  write ( *, '(2x,a4,4x,g14.6)' ) name, result
 
  name = 'X**2'
  call polygon_xx_2d ( npts, x, y, result )
  write ( *, '(2x,a4,4x,g14.6)' ) name, result
 
  name = 'XY'
  call polygon_xy_2d ( npts, x, y, result )
  write ( *, '(2x,a4,4x,g14.6)' ) name, result
 
  name = 'Y**2'
  call polygon_yy_2d ( npts, x, y, result )
  write ( *, '(2x,a4,4x,g14.6)' ) name, result
 
  return
end
subroutine test25

!*****************************************************************************80
!
!! TEST25 tests PYRAMID_UNIT_O**_3D, PYRAMID_VOLUME_3D.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: rule_num = 10

  character ( len = 7 ), allocatable, dimension ( : ) :: column_label
  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) jhi
  integer   ( kind = 4 ) jlo
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) pyramid_unit_volume_3d
  real      ( kind = 8 ), allocatable, dimension ( :, : ) :: result
  integer   ( kind = 4 ), dimension ( rule_num ) :: order = (/ &
    1, 5, 6, 8, 8, 9, 13, 18, 27, 48 /)

  num = function_3d_num ( )

  allocate ( column_label(1:num) )
  allocate ( result(1:rule_num,1:num) )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST25'
  write ( *, '(a)' ) '  For the unit pyramid, we approximate integrals with:'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O01_3D, a 1 point rule.'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O05_3D, a 5 point rule.'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O06_3D, a 6 point rule.'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O08_3D, an 8 point rule.'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O08b_3D, an 8 point rule.'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O09_3D, a 9 point rule.'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O13_3D, a 13 point rule.'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O18_3D, a 18 point rule.'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O27_3D, a 27 point rule.'
  write ( *, '(a)' ) '  PYRAMID_UNIT_O48_3D, a 48 point rule.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  PYRAMID_UNIT_VOLUME_3D computes the volume of a unit pyramid.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Volume = ', pyramid_unit_volume_3d ( )
  write ( *, '(a)' ) ' '
    
  do j = 1, num

    call function_3d_set ( 'SET', j )

    column_label(j) = function_3d_name(j)

    call pyramid_unit_o01_3d ( function_3d, result(1,j) )
    call pyramid_unit_o05_3d ( function_3d, result(2,j) )
    call pyramid_unit_o06_3d ( function_3d, result(3,j) )
    call pyramid_unit_o08_3d ( function_3d, result(4,j) )
    call pyramid_unit_o08b_3d ( function_3d, result(5,j) )
    call pyramid_unit_o09_3d ( function_3d, result(6,j) )
    call pyramid_unit_o13_3d ( function_3d, result(7,j) )
    call pyramid_unit_o18_3d ( function_3d, result(8,j) )
    call pyramid_unit_o27_3d ( function_3d, result(9,j) )
    call pyramid_unit_o48_3d ( function_3d, result(10,j) )

  end do

  do jlo = 1, num, 5
    jhi = min ( jlo + 4, num )
    write ( *, '(a)' ) ' '
    write ( *, '(2x,a5,3x,a7,7x,a7,7x,a7,7x,a7,7x,a7)' ) &
      'Order', column_label(jlo:jhi)
    write ( *, '(a)' ) ' '
    do i = 1, rule_num
      write ( *, '(2x,i5,5g14.6)' ) order(i), result(i,jlo:jhi)
    end do
  end do

  deallocate ( column_label)
  deallocate ( result )

  return
end
subroutine test255

!*****************************************************************************80
!
!! TEST255 tests PYRAMID_UNIT_MONOMIAL_3D.
!
!  Modified:
!
!    24 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ) alpha
  integer   ( kind = 4 ) beta
  integer   ( kind = 4 ) :: degree_max = 4
  integer   ( kind = 4 ) gamma
  real      ( kind = 8 ) pyramid_unit_monomial_3d
  real      ( kind = 8 ) pyramid_unit_volume_3d
  real      ( kind = 8 ) value

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST255'
  write ( *, '(a)' ) '  For the unit pyramid,'
  write ( *, '(a)' ) '  PYRAMID_UNIT_MONOMIAL_3D returns the exact value of the'
  write ( *, '(a)' ) ' integral of X^ALPHA Y^BETA Z^GAMMA'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Volume = ', pyramid_unit_volume_3d ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     ALPHA      BETA     GAMMA      INTEGRAL'
  write ( *, '(a)' ) ' '

  do alpha = 0, degree_max
    do beta = 0, degree_max - alpha
      do gamma = 0, degree_max - alpha - beta
        value = pyramid_unit_monomial_3d ( alpha, beta, gamma )
        write ( *, '(2x,i8,2x,i8,2x,i8,2x,g14.6)' ) alpha, beta, gamma, value
      end do
    end do
  end do

  return
end
subroutine test26

!*****************************************************************************80
!
!! TEST26 tests QMULT_1D.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ) a
  real      ( kind = 8 ) b
  character ( len = 7 ) function_1d_name
  real      ( kind = 8 ), external :: function_1d
  integer   ( kind = 4 ) function_1d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) qmult_1d
  real      ( kind = 8 ) result

  a = -1.0D+00
  b = 1.0D+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST26'
  write ( *, '(a)' ) '  QMULT_1D approximates an integral on a'
  write ( *, '(a)' ) '    one-dimensional interval.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We use the interval:'
  write ( *, '(a,g14.6)' ) '  A = ', a
  write ( *, '(a,g14.6)' ) '  B = ', b
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    F(X)     QMULT_1D'
  write ( *, '(a)' ) ' '
 
  num = function_1d_num ( )

  do i = 1, num

    call function_1d_set ( 'SET', i )

    result = qmult_1d ( function_1d, a, b )
    write ( *, '(2x,a7,f14.8)' ) function_1d_name ( i ), result
 
  end do

  return
end
subroutine test27

!*****************************************************************************80
!
!! TEST27 tests SIMPLEX_ND.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) function_nd_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result
  real      ( kind = 8 ), allocatable, dimension ( :, : ) :: v

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST27'
  write ( *, '(a)' ) '  SIMPLEX_ND approximates integrals inside an' 
  write ( *, '(a)' ) '    arbitrary simplex in ND.'
  write ( *, '(a)' ) ' '
 
  do n = 2, 4
 
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Spatial dimension N = ', n

    allocate ( v(1:n,1:n+1) )
!
!  Restore values of simplex.
!
    call setsim ( n, v )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Simplex vertices:'
    write ( *, '(a)' ) ' '
 
    do i = 1, n
      write ( *, '(5f4.0)' ) v(i,1:n+1)
    end do
 
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  F(X)    SIMPLEX_ND'
    write ( *, '(a)' ) ' '

    num = function_nd_num ( )

    do i = 1, num

      call function_nd_set ( 'SET', i )

      call simplex_nd ( function_nd, n, v, result )
      write ( *, '(2x,a7,f14.8)' ) function_nd_name(i), result
 
      call setsim ( n, v )

    end do

    deallocate ( v )

  end do
 
  return
end
subroutine test28

!*****************************************************************************80
!
!! TEST28 tests SIMPLEX_VOLUME_ND.
!
!  Modified:
!
!    06 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) i
  integer ( kind = 4 ) n
  real    ( kind = 8 ) simplex_volume_nd
  real    ( kind = 8 ), allocatable, dimension ( :, : ) :: v
  real    ( kind = 8 ) volume

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST28'
  write ( *, '(a)' ) '  SIMPLEX_VOLUME_ND computes the volume of a simplex'
  write ( *, '(a)' ) '    in N dimensions.'
  write ( *, '(a)' ) ' '

  do n = 2, 4

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Spatial dimension N = ', n

    allocate ( v(1:n,1:n+1) )
!
!  Set the values of the simplex.
!
    call setsim ( n, v )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Simplex vertices:'
    write ( *, '(a)' ) ' '

    do i = 1, n
      write ( *, '(5f4.0)' ) v(i,1:n+1)
    end do

    volume = simplex_volume_nd ( n, v )

    write ( *, '(a)' ) ' '
    write ( *, '(a,g14.6)' ) '  Volume is ', volume

    deallocate ( v )

  end do

  return
end
subroutine test29

!*****************************************************************************80
!
!! TEST29 tests SIMPLEX_UNIT_**_ND.
!
!  Modified:
!
!    10 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) n
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) result3
  real      ( kind = 8 ) result4
  real      ( kind = 8 ) simplex_unit_volume_nd
  real      ( kind = 8 ) volume

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST29'
  write ( *, '(a)' ) '  For integrals in the unit simplex in ND,'
  write ( *, '(a)' ) '  SIMPLEX_UNIT_01_ND uses a formula of degree 1.'
  write ( *, '(a)' ) '  SIMPLEX_UNIT_03_ND uses a formula of degree 3.'
  write ( *, '(a)' ) '  SIMPLEX_UNIT_05_ND uses a formula of degree 5.'
  write ( *, '(a)' ) '  SIMPLEX_UNIT_05_2_ND uses a formula of degree 5.'

  do i = 1, 6

    call function_nd_set ( 'SET', i )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Check the integral of ', function_nd_name ( i )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  N     Volume         #1              #3' // &
      '              #5              #5.2'
    write ( *, '(a)' ) ' '

    do n = 2, 16

      call simplex_unit_01_nd ( function_nd, n, result1 )
      call simplex_unit_03_nd ( function_nd, n, result2 )
      call simplex_unit_05_nd ( function_nd, n, result3 )
      call simplex_unit_05_2_nd ( function_nd, n, result4 )

      volume = simplex_unit_volume_nd ( n )
      write ( *, '(2x,i2,2x,g13.5,2x,g13.5,2x,g13.5,2x,g13.5,2x,g13.5)' ) &
        n, volume, result1, result2, result3, result4

    end do

  end do

  return
end
subroutine test30

!*****************************************************************************80
!
!! TEST30 tests SPHERE_UNIT_**_3D.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) result3
  real      ( kind = 8 ) result4
  real      ( kind = 8 ) sphere_unit_area_nd

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST30'
  write ( *, '(a)' ) '  For integrals on the unit sphere in 3D:'
  write ( *, '(a)' ) '  SPHERE_UNIT_07_3D uses a formula of degree 7.'
  write ( *, '(a)' ) '  SPHERE_UNIT_11_3D uses a formula of degree 11.'
  write ( *, '(a)' ) '  SPHERE_UNIT_14_3D uses a formula of degree 14.'
  write ( *, '(a)' ) '  SPHERE_UNIT_15_3D uses a formula of degree 15.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Unit sphere area = ', sphere_unit_area_nd ( 3 )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    F(X)    S3S07        S3S11         S3S14' // &
    '         S3S15      '
  write ( *, '(a)' ) ' '
 
  num = function_3d_num ( )

  do i = 1, num

    call function_3d_set ( 'SET', i )

    call sphere_unit_07_3d ( function_3d, result1 ) 
    call sphere_unit_11_3d ( function_3d, result2 )
    call sphere_unit_14_3d ( function_3d, result3 )
    call sphere_unit_15_3d ( function_3d, result4 )
 
    write ( *, '(2x,a7,4f14.8)' ) function_3d_name(i), result1, result2, result3, result4

  end do
 
  return
end
subroutine test31

!*****************************************************************************80
!
!! TEST31 tests SPHERE_UNIT_**_ND.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) function_nd_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) result3
  real      ( kind = 8 ) result4
  real      ( kind = 8 ) result5
  real      ( kind = 8 ) result6
  real      ( kind = 8 ) sphere_unit_area_nd

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST31'
  write ( *, '(a)' ) '  For integrals on the unit sphere in ND:'
  write ( *, '(a)' ) '  SPHERE_UNIT_03_ND uses a formula of degree 3;'
  write ( *, '(a)' ) '  SPHERE_UNIT_04_ND uses a formula of degree 4;'
  write ( *, '(a)' ) '  SPHERE_UNIT_05_ND uses a formula of degree 5.'
  write ( *, '(a)' ) '  SPHERE_UNIT_07_1_ND uses a formula of degree 7.'
  write ( *, '(a)' ) '  SPHERE_UNIT_07_2_ND uses a formula of degree 7.'
  write ( *, '(a)' ) '  SPHERE_UNIT_11_ND uses a formula of degree 11.'
  write ( *, '(a)' ) ' '

  do n = 3, 10

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Spatial dimension N = ', n
    write ( *, '(a,g14.6)' ) '  Unit sphere area = ', sphere_unit_area_nd ( n )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '    Rule:     #3            #4            #5'          
    write ( *, '(a)' ) '              #7.1          #7.2          #11'
    write ( *, '(a)' ) '    Function'
    write ( *, '(a)' ) ' '

    num = function_nd_num ( )

    do i = 1, num

      call function_nd_set ( 'SET', i )

      call sphere_unit_03_nd ( function_nd, n, result1 )
      call sphere_unit_04_nd ( function_nd, n, result2 )
      call sphere_unit_05_nd ( function_nd, n, result3 )
      call sphere_unit_07_1_nd ( function_nd, n, result4 )
      call sphere_unit_07_2_nd ( function_nd, n, result5 )
      call sphere_unit_11_nd ( function_nd, n, result6 )

      write ( *, '(2x,a7,3f14.8)' ) function_nd_name(i), result1, result2, result3
      write ( *, '(2x,7x,3f14.8)' )           result4, result5, result6

    end do

  end do

  return
end
subroutine test32

!*****************************************************************************80
!
!! TEST32 tests SPHERE_05_ND, SPHERE_07_1_ND.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: n_max = 5

  real      ( kind = 8 ) center(n_max)
  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) function_nd_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) r
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) sphere_area_nd

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST32'
  write ( *, '(a)' ) '  For integrals on a sphere in ND:'
  write ( *, '(a)' ) '  SPHERE_05_ND uses a formula of degree 5.'
  write ( *, '(a)' ) '  SPHERE_07_1_ND uses a formula of degree 7.'
  write ( *, '(a)' ) ' '
  r = 2.0D+00
  center(1:n_max) = 1.0D+00

  do n = 2, 4

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Spatial dimension N = ', n
    write ( *, '(a)' ) '  Sphere center = '
    write ( *, '(5g14.6)' ) center(1:n)
    write ( *, '(a,g14.6)' ) '  Sphere radius = ', r
    write ( *, '(a,g14.6)' ) '  Sphere area = ', sphere_area_nd ( n, r )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '    Rule:     #5           #7.1'
    write ( *, '(a)' ) '    Function'
    write ( *, '(a)' ) ' '

    num = function_nd_num ( )

    do i = 1, num

      call function_nd_set ( 'SET', i )

      call sphere_05_nd ( function_nd, n, center, r, result1 )
      call sphere_07_1_nd ( function_nd, n, center, r, result2 )

      write ( *, '(2x,a7,2f14.8)' ) function_nd_name(i), result1, result2

    end do

  end do

  return
end
subroutine test322

!*****************************************************************************80
!
!! TEST322 tests SPHERE_CAP_AREA_3D, SPHERE_CAP_AREA_ND.
!
!  Modified:
!
!    17 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: dim_num = 3

  real    ( kind = 8 ) area1
  real    ( kind = 8 ) area2
  real    ( kind = 8 ) center(dim_num)
  real    ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: test_num = 12
  real    ( kind = 8 ) :: r = 1.0D+00
  real    ( kind =8 ) sphere_area_3d

  center(1:3) = (/ 0.0D+00, 0.0D+00, 0.0D+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST322'
  write ( *, '(a)' ) '  SPHERE_CAP_AREA_3D computes the volume of a'
  write ( *, '(a)' ) '    3D spherical cap, defined by a plane that cuts the'
  write ( *, '(a)' ) '    sphere to a thickness of H units.'
  write ( *, '(a)' ) '  SPHERE_CAP_AREA_ND computes the volume of an'
  write ( *, '(a)' ) '    ND spherical cap, defined by a plane that cuts the'
  write ( *, '(a)' ) '    sphere to a thickness of H units.'

  area1 = sphere_area_3d ( r )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Area of the total sphere in 3D = ', area1

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '        R           H           Cap         Cap'
  write ( *, '(a)' ) '                                area_3d     area_nd'
  write ( *, '(a)' ) ' '

  do i = 0, test_num+1

    h = 2.0D+00 * r * real ( i, kind = 8 ) / real ( test_num, kind = 8 )

    call sphere_cap_area_3d ( r, h, area1 )

    call sphere_cap_area_nd ( dim_num, r, h, area2 )

    write ( *, '(2x,5f12.6)' ) r, h, area1, area2

  end do

  return
end
subroutine test324

!*****************************************************************************80
!
!! TEST324 tests SPHERE_CAP_VOLUME_2D, SPHERE_CAP_VOLUME_ND.
!
!  Modified:
!
!    17 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: dim_num = 2

  real    ( kind = 8 ) center(dim_num)
  real    ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: test_num = 12
  real    ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real    ( kind = 8 ) :: r = 1.0D+00
  real    ( kind = 8 ) volume1
  real    ( kind = 8 ) volume2

  center(1:2) = (/ 0.0D+00, 0.0D+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST324'
  write ( *, '(a)' ) '  SPHERE_CAP_VOLUME_2D computes the volume (area) of a'
  write ( *, '(a)' ) '    spherical cap, defined by a plane that cuts the'
  write ( *, '(a)' ) '    sphere to a thickness of H units.'
  write ( *, '(a)' ) '  SPHERE_CAP_VOLUME_ND does the same operation,'
  write ( *, '(a)' ) '    but in N dimensions.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Using a radius R = ', r

  call sphere_volume_2d ( r, volume1 )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Volume of the total sphere in 2D = ', volume1

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '        H           Cap        Cap'
  write ( *, '(a)' ) '                    vol_2d     vol_nd'
  write ( *, '(a)' ) ' '

  do i = 0, test_num+1

    h = 2.0D+00 * r * real ( i, kind = 8 ) / real ( test_num, kind = 8 )

    call sphere_cap_volume_2d ( r, h, volume1 )

    call sphere_cap_volume_nd ( dim_num, r, h, volume2 )

    write ( *, '(2x,3f12.6)' ) h, volume1, volume2

  end do

  return
end
subroutine test326

!*****************************************************************************80
!
!! TEST326 tests SPHERE_CAP_VOLUME_3D, SPHERE_CAP_VOLUME_ND.
!
!  Modified:
!
!    17 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ), parameter :: dim_num = 3

  real    ( kind = 8 ) center(dim_num)
  real    ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: test_num = 12
  real    ( kind = 8 ) :: r = 1.0D+00
  real    ( kind = 8 ) volume1
  real    ( kind = 8 ) volume2

  center(1:3) = (/ 0.0D+00, 0.0D+00, 0.0D+00 /)
  
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST326'
  write ( *, '(a)' ) '  SPHERE_CAP_VOLUME_3D computes the volume of a'
  write ( *, '(a)' ) '    spherical cap, defined by a plane that cuts the'
  write ( *, '(a)' ) '    sphere to a thickness of H units.'
  write ( *, '(a)' ) '  SPHERE_CAP_VOLUME_ND does the same operation,'
  write ( *, '(a)' ) '    but in N dimensions.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Using a radius R = ', r

  call sphere_volume_3d ( r, volume1 )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Volume of the total sphere in 3D = ', volume1

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '        H           Cap        Cap'
  write ( *, '(a)' ) '                    volume_3d  volume_nd'
  write ( *, '(a)' ) ' '

  do i = 0, test_num+1

    h = 2.0D+00 * r * real ( i, kind = 8 ) / real ( test_num, kind = 8 )

    call sphere_cap_volume_3d ( r, h, volume1 )

    call sphere_cap_volume_nd ( dim_num, r, h, volume2 )

    write ( *, '(2x,3f12.6)' ) h, volume1, volume2

  end do

  return
end
subroutine test33

!*****************************************************************************80
!
!! TEST33 tests SPHERE_CAP_AREA_ND, SPHERE_CAP_VOLUME_ND.
!
!  Modified:
!
!    17 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real    ( kind = 8 ) area
  real    ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: n = 12
  integer ( kind = 4 ) dim_num
  real    ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real    ( kind = 8 ) r
  real    ( kind = 8 ) sphere_area_nd
  real    ( kind = 8 ) volume

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST33'
  write ( *, '(a)' ) '  For a sphere in ND:'
  write ( *, '(a)' ) '  SPHERE_CAP_AREA_ND computes the area '
  write ( *, '(a)' ) '    of a spherical cap.'
  write ( *, '(a)' ) '  SPHERE_CAP_VOLUME_ND computes the volume '
  write ( *, '(a)' ) '    of a spherical cap.'
  write ( *, '(a)' ) ' '

  r = 1.0D+00

  do dim_num = 2, 5

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Spatial dimension N = ', dim_num
    write ( *, '(a,g14.6)' ) '  Radius =       ', r
    write ( *, '(a,g14.6)' ) '  Area =         ', sphere_area_nd ( dim_num, r )
    call sphere_volume_nd ( dim_num, r, volume )
    write ( *, '(a,g14.6)' ) '  Volume =       ', volume

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '                 Sphere         Sphere' 
    write ( *, '(a)' ) '                 cap            cap'
    write ( *, '(a)' ) '      H          area           volume'
    write ( *, '(a)' ) ' '

    do i = 0, n + 1
      h = real ( 2 * i, kind = 8 ) * r / real ( n, kind = 8 )
      call sphere_cap_area_nd ( dim_num, r, h, area )
      call sphere_cap_volume_nd ( dim_num, r, h, volume )
      write ( *, '(2x,f8.4,2x,g14.6,2x,g14.6)' ) h, area, volume
    end do

  end do

  return
end
subroutine test335

!*****************************************************************************80
!
!! TEST335 tests SPHERE_SHELL_03_ND.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: n_max = 3

  real      ( kind = 8 ) center(n_max)
  real      ( kind = 8 ), external :: function_nd
  character ( len = 7 ) function_nd_name
  integer   ( kind = 4 ) function_nd_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result3
  real      ( kind = 8 ) result4
  real      ( kind = 8 ) result5
  real      ( kind = 8 ) result6
  real      ( kind = 8 ) r1
  real      ( kind = 8 ) r2
  real      ( kind = 8 ) sphere_shell_volume_nd

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST335'
  write ( *, '(a)' ) '  For integrals inside a spherical shell in ND:'
  write ( *, '(a)' ) '  SPHERE_SHELL_03_ND approximates the integral.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We compare these results with those computed by'
  write ( *, '(a)' ) '  from the difference of two ball integrals:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  BALL_F1_ND approximates the integral;'
  write ( *, '(a)' ) '  BALL_F3_ND approximates the integral.'
  write ( *, '(a)' ) ' '

  do j = 1, 2

    if ( j == 1 ) then
      r1 = 0.0D+00
      r2 = 1.0D+00
      center(1:n_max) = 0.0D+00
    else
      r1 = 2.0D+00
      r2 = 3.0D+00
      center(1:n_max) = (/ 1.0D+00, -1.0D+00, 2.0D+00 /)
    end if

    do n = 2, n_max

      write ( *, '(a)' ) ' '
      write ( *, '(a,i8)' ) '  Spatial dimension N = ', n
      write ( *, '(a)' ) '  Sphere center:'
      write ( *, '(2x,3f10.4)' ) center(1:n)
      write ( *, '(a,g14.6)' ) '  Inner sphere radius = ', r1
      write ( *, '(a,g14.6)' ) '  Outer sphere radius = ', r2
      write ( *, '(a,g14.6)' ) '  Spherical shell volume = ', &
        sphere_shell_volume_nd ( n, r1, r2 )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '    Rule:      #3       F1(R2)-F1(R1)  ' // &      
        'F3(R2)-F3(R1)'
      write ( *, '(a)' ) '    F(X)'
      write ( *, '(a)' ) ' '

      num = function_nd_num ( )

      do i = 1, num

        call function_nd_set ( 'SET', i )

        call sphere_shell_03_nd ( function_nd, n, center, r1, r2, result1 )

        call ball_f1_nd ( function_nd, n, center, r1, result3 )
        call ball_f1_nd ( function_nd, n, center, r2, result4 )
      
        call ball_f3_nd ( function_nd, n, center, r1, result5 )
        call ball_f3_nd ( function_nd, n, center, r2, result6 )

        write ( *, '(2x,a7,4f14.8)' ) function_nd_name(i), result1, &
          result4-result3, result6-result5

      end do

    end do

  end do

  return
end
subroutine test34

!*****************************************************************************80
!
!! TEST34 tests SPHERE_UNIT_AREA_ND, SPHERE_UNIT_AREA_VALUES.
!
!  Modified:
!
!    17 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real    ( kind = 8 ) area
  real    ( kind = 8 ) area2
  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) n_data
  real    ( kind = 8 ) sphere_unit_area_nd

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST34:'
  write ( *, '(a)' ) '  SPHERE_UNIT_AREA_ND evaluates the area of the unit'
  write ( *, '(a)' ) '  sphere in N dimensions.'
  write ( *, '(a)' ) '  SPHERE_UNIT_AREA_VALUES returns some test values.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     dim_num    Exact          Computed'
  write ( *, '(a)' ) '             Area           Area'
  write ( *, '(a)' ) ' '

  n_data = 0

  do

    call sphere_unit_area_values ( n_data, dim_num, area )

    if ( n_data == 0 ) then
      exit
    end if

    area2 = sphere_unit_area_nd ( dim_num )

    write ( *, '(2x,i8,2x,f10.6,2x,f10.6)' ) dim_num, area, area2

  end do

  return
end
subroutine test345

!*****************************************************************************80
!
!! TEST345 tests SPHERE_UNIT_VOLUME_ND, SPHERE_UNIT_VOLUME_VALUES.
!
!  Modified:
!
!    17 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) n_data
  real    ( kind = 8 ) sphere_unit_volume_nd
  real    ( kind = 8 ) volume
  real    ( kind = 8 ) volume2

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST345:'
  write ( *, '(a)' ) '  SPHERE_UNIT_VOLUME_ND evaluates the area of the unit'
  write ( *, '(a)' ) '  sphere in N dimensions.'
  write ( *, '(a)' ) '  SPHERE_UNIT_VOLUME_VALUES returns some test values.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     dim_num    Exact          Computed'
  write ( *, '(a)' ) '             Volume         Volume'
  write ( *, '(a)' ) ' '

  n_data = 0

  do

    call sphere_unit_volume_values ( n_data, dim_num, volume )

    if ( n_data == 0 ) then
      exit
    end if

    volume2 = sphere_unit_volume_nd ( dim_num )

    write ( *, '(2x,i8,2x,f10.6,2x,f10.6)' ) dim_num, volume, volume2

  end do

  return
end
subroutine test35

!*****************************************************************************80
!
!! TEST35 tests SQUARE_UNIT_SET, RECTANGLE_SUB_2D.
!
!  Modified:
!
!    15 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  integer   ( kind = 4 ) nsub(2)
  real      ( kind = 8 ) result
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ) xval(2)
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab
  real      ( kind = 8 ) yval(2)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST35'
  write ( *, '(a)' ) '  SQUARE_UNIT_SET sets up a quadrature rule '
  write ( *, '(a)' ) '    on a unit square.'
  write ( *, '(a)' ) '  RECTANGLE_SUB_2D applies it to subrectangles of an'
  write ( *, '(a)' ) '    arbitrary rectangle.'
  write ( *, '(a)' ) ' '
!
!  Set the location of the square.
!
  xval(1) = 1.0D+00
  yval(1) = 2.0D+00

  xval(2) = 3.0D+00
  yval(2) = 3.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The corners of the rectangle are:'
  write ( *, '(a)' ) ' '
  write ( *, '(2g14.6)' ) xval(1), yval(1)
  write ( *, '(2g14.6)' ) xval(2), yval(2)
!
!  Get the quadrature abscissas and weights for a unit square.
!
  rule = 2
  call square_unit_size ( rule, order )

  allocate ( xtab(1:order) )
  allocate ( ytab(1:order) )
  allocate ( weight(1:order) )

  call square_unit_set ( rule, order, xtab, ytab, weight )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Using unit square integration rule number ', rule
  write ( *, '(a,i8)' ) '  Order of rule is ', order
!
!  Set the function.
!
  num = function_2d_num ( )

  do i = 1, num

    call function_2d_set ( 'SET', i )
!
!  Try an increasing number of subdivisions.
!
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '    Function  Subdivisions  Integral'
    write ( *, '(a)' ) ' '

    do j = 1, 5

      nsub(1) = j
      nsub(2) = 2 * j

      call rectangle_sub_2d ( function_2d, xval, yval, nsub, order, xtab, ytab, &
        weight, result )

      write ( *, '(2x,a7,2i4, f14.8)' ) function_2d_name(i), nsub(1), nsub(2), result

    end do

  end do

  deallocate ( weight )
  deallocate ( xtab )
  deallocate ( ytab )

  return
end
subroutine test36

!*****************************************************************************80
!
!! TEST36 tests SQUARE_UNIT_SET and SQUARE_SUM.
!
!  Modified:
!
!    15 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: dim_num = 2
  integer   ( kind = 4 ), parameter :: rule_max = 6

  real      ( kind = 8 ), dimension ( dim_num ) :: center = (/ 2.0D+00, 2.0D+00 /)
  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) r
  real      ( kind = 8 ) result(rule_max)
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab

  r = 3.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST36'
  write ( *, '(a)' ) '  SQUARE_UNIT_SET sets up quadrature on the unit square;'
  write ( *, '(a)' ) '  SQUARE_SUM carries it out on an arbitrary square.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Square center:'
  write ( *, '(a,2g14.6)' ) '  CENTER = ', center(1:dim_num)
  write ( *, '(a,g14.6)' ) '  Square radius is ', r

  do ilo = 1, rule_max, 5

    ihi = min ( ilo + 4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a,5(i6,7x))' ) 'Rule:  ', ( rule, rule = ilo, ihi )
    write ( *, '(a)' ) '  Function '
    write ( *, '(a)' ) ' '

    num = function_2d_num ( )

    do i = 1, num

      call function_2d_set ( 'SET', i )

      do rule = ilo, ihi

        call square_unit_size ( rule, order )

        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( weight(1:order) )

        call square_unit_set ( rule, order, xtab, ytab, weight )

        call square_sum ( function_2d, center, r, order, xtab, ytab, weight, &
          result(rule) )

        deallocate ( weight )
        deallocate ( xtab )
        deallocate ( ytab )

      end do

      write ( *, '(2x,a7,5f13.7)' ) function_2d_name(i), result(ilo:ihi)

    end do

  end do
 
  return
end
subroutine test37

!*****************************************************************************80
!
!! TEST37 tests SQUARE_UNIT_SET and SQUARE_UNIT_SUM.
!
!  Modified:
!
!    15 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: rule_max = 6

  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) result(rule_max)
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST37'
  write ( *, '(a)' ) '  SQUARE_UNIT_SET sets up quadrature on the unit square;'
  write ( *, '(a)' ) '  SQUARE_UNIT_SUM carries it out on the unit square.'
  write ( *, '(a)' ) ' '
 
  do ilo = 1, rule_max, 5

    ihi = min ( ilo + 4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a,5i6)' ) 'Rule:    ', ( rule, rule = ilo, ihi )
    write ( *, '(2x,a)' ) 'Function '
    write ( *, '(a)' ) ' '

    num = function_2d_num ( )

    do i = 1, num

      call function_2d_set ( 'SET', i )

      do rule = ilo, ihi

        call square_unit_size ( rule, order )

        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( weight(1:order) )

        call square_unit_set ( rule, order, xtab, ytab, weight )

        call square_unit_sum ( function_2d, order, xtab, ytab, weight, &
          result(rule) )

        deallocate ( weight )
        deallocate ( xtab )
        deallocate ( ytab )

      end do

      write ( *, '(2x,a7,5f13.7)' ) function_2d_name(i), result(ilo:ihi)

    end do

  end do
 
  return
end
subroutine test38

!*****************************************************************************80
!
!! TEST38 tests TETRA_07, TETRA_TPRODUCT.
!
!  Modified:
!
!    05 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: order_max = 9

  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) result3(order_max)
  real      ( kind = 8 ) tetra_unit_volume
  real      ( kind = 8 ) tetra_volume
  real      ( kind = 8 ), dimension ( 4 ) :: x = (/ &
    1.0D+00, 4.0D+00, 1.0D+00, 1.0D+00 /)
  real      ( kind = 8 ), dimension ( 4 ) :: y = (/ &
    2.0D+00, 2.0D+00, 3.0D+00, 2.0D+00 /)
  real      ( kind = 8 ), dimension ( 4 ) :: z = (/ &
    6.0D+00, 6.0D+00, 6.0D+00, 8.0D+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST38'
  write ( *, '(a)' ) '  For integrals inside an arbitrary tetrahedron:'
  write ( *, '(a)' ) '  TETRA_07 uses a formula of degree 7;'
  write ( *, '(a)' ) '  TETRA_TPRODUCT uses a triangular product formula '
  write ( *, '(a)' ) '    of varying degree.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Tetrahedron vertices:'
  write ( *, '(a)' ) ' '
  do i = 1, 4
    write ( *, '(2x,3f4.0)' ) x(i), y(i), z(i)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Tetrahedron unit volume = ', tetra_unit_volume ( )
  write ( *, '(a,g14.6)' ) '  Tetrahedron Volume = ', tetra_volume ( x, y, z )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  F(X)    TETRA_07'
  write ( *, '(a)' ) '          TETRA_TPRODUCT(1:4)'
  write ( *, '(a)' ) '          TETRA_TPRODUCT(5:8)'
  write ( *, '(a)' ) '          TETRA_TPRODUCT(9)'
  write ( *, '(a)' ) ' '

  num = function_3d_num ( )

  do i = 1, num

    call function_3d_set ( 'SET', i )

    call tetra_07 ( function_3d, x, y, z, result2 )

    do order = 1, order_max
      call tetra_tproduct ( function_3d, order, x, y, z, result3(order) )
    end do

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a7,4f16.10)' ) function_3d_name(i), result2
    write ( *, '(2x,7x,4f16.10)' )                  result3(1:4)
    write ( *, '(2x,7x,4f16.10)' )                  result3(5:8)
    write ( *, '(2x,7x,4f16.10)' )                  result3(9)

  end do

  return
end
subroutine test39

!*****************************************************************************80
!
!! TEST39 tests TETRA_UNIT_SET and TETRA_UNIT_SUM.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: rule_max = 8

  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) result(rule_max)
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ztab

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST39'
  write ( *, '(a)' ) '  TETRA_UNIT_SET sets quadrature rules '
  write ( *, '(a)' ) '    for the unit tetrahedron;'
  write ( *, '(a)' ) '  TETRA_UNIT_SUM applies them to the unit tetrahedron.'
  write ( *, '(a)' ) ' '

  do ilo = 1, rule_max, 5

    ihi = min ( ilo +  4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a,5i6)' ) 'Rule:   ', ( rule, rule = ilo, ihi )
    write ( *, '(2x,a)' ) 'Function'
    write ( *, '(a)' ) ' '

    num = function_3d_num ( )

    do i = 1, num

      call function_3d_set ( 'SET', i )

      do rule = ilo, ihi

        call tetra_unit_size ( rule, order )

        allocate ( weight(1:order) )
        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( ztab(1:order) )

        call tetra_unit_set ( rule, order, xtab, ytab, ztab, weight )
 
        call tetra_unit_sum ( function_3d, order, xtab, ytab, ztab, weight, &
          result(rule) )

        deallocate ( weight )
        deallocate ( xtab )
        deallocate ( ytab )
        deallocate ( ztab )

      end do

      write ( *, '(2x,a7,5f14.6)' ) function_3d_name(i), result(ilo:ihi)

    end do

  end do

  return
end
subroutine test40

!*****************************************************************************80
!
!! TEST40 tests TETRA_UNIT_SET and TETRA_SUM.
!
!  Modified:
!
!    02 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: rule_max = 8

  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) result(rule_max)
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ) value
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), dimension ( 4 ) :: x = (/ &
    1.0D+00, 4.0D+00, 1.0D+00, 1.0D+00 /)
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ), dimension ( 4 ) :: y = (/ &
    2.0D+00, 2.0D+00, 3.0D+00, 2.0D+00 /)
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab
  real      ( kind = 8 ), dimension ( 4 ) :: z = (/ &
    6.0D+00, 6.0D+00, 6.0D+00, 8.0D+00 /)
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ztab

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST40'
  write ( *, '(a)' ) '  TETRA_UNIT_SET sets quadrature rules '
  write ( *, '(a)' ) '    for the unit tetrahedron;'
  write ( *, '(a)' ) '  TETRA_SUM applies them to an arbitrary tetrahedron.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Tetrahedron vertices:'
  write ( *, '(a)' ) ' '
  do i = 1, 4
    write ( *, '(2x,3f6.2)' ) x(i), y(i), z(i)
  end do

  do ilo = 1, rule_max, 5

    ihi = min ( ilo +  4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a,5(i7,7x))' ) 'Rule:    ', ( rule, rule = ilo, ihi )
    write ( *, '(2x,a)' ) 'Function '
    write ( *, '(a)' ) ' '
!
!  Believe it or not, if I remove these unused debugging statements,
!  the program will fail under the G95 compiler.  Somehow, somewhere,
!  there must be a memory error, but I have not figured out where it is.
!
    num = function_3d_num ( )

    do i = 1, num

      call function_3d_set ( 'SET', i )

      do rule = ilo, ihi

        call tetra_unit_size ( rule, order )

        allocate ( weight(1:order) )
        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( ztab(1:order) )

        call tetra_unit_set ( rule, order, xtab, ytab, ztab, weight )

        call tetra_sum ( function_3d, x, y, z, order, xtab, ytab, ztab, weight, &
          value )

        result(rule) = value

        deallocate ( weight )
        deallocate ( xtab )
        deallocate ( ytab )
        deallocate ( ztab )

      end do

      write ( *, '(2x,a7,5g14.7)' ) function_3d_name(i), result(ilo:ihi)

    end do

  end do

  return
end
subroutine test41

!*****************************************************************************80
!
!! TEST41 tests TRIANGLE_UNIT_SET, TRIANGLE_SUB.
!
!  Discussion:
!
!    Break up the triangle into NSUB*NSUB equal subtriangles.  Approximate 
!    the integral over the triangle by the sum of the integrals over each
!    subtriangle.
!
!  Modified:
!
!    07 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  integer   ( kind = 4 ) nsub
  real      ( kind = 8 ) result
  integer   ( kind = 4 ) rule
  integer   ( kind = 4 ) triangle_unit_size
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ), dimension ( 3 ) :: xval = (/ &
    0.0D+00, 0.0D+00, 1.0D+00 /)
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab
  real      ( kind = 8 ), dimension ( 3 ) :: yval = (/ &
    0.0D+00, 1.0D+00, 0.0D+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST41'
  write ( *, '(a)' ) '  TRIANGLE_UNIT_SET sets up a quadrature rule '
  write ( *, '(a)' ) '    on a triangle.'
  write ( *, '(a)' ) '  TRIANGLE_SUB applies it to subtriangles of an'
  write ( *, '(a)' ) '    arbitrary triangle.'
  write ( *, '(a)' ) ' '

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Triangle vertices:'
  write ( *, '(a)' ) ' '
  write ( *, '(2x,2g14.6)' ) xval(1), yval(1)
  write ( *, '(2x,2g14.6)' ) xval(2), yval(2)
  write ( *, '(2x,2g14.6)' ) xval(3), yval(3)
!
!  Get the quadrature abscissas and weights for a unit triangle.
!
  rule = 3
  order = triangle_unit_size ( rule )

  allocate ( xtab(1:order) )
  allocate ( ytab(1:order) )
  allocate ( weight(1:order) )

  call triangle_unit_set ( rule, order, xtab, ytab, weight )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Using unit triangle quadrature rule ', rule
  write ( *, '(a,i8)' ) '  Rule order = ', order
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Function Nsub  Result'
  write ( *, '(a)' ) ' '
!
!  Set the function.
!
  num = function_2d_num ( )

  do i = 1, num

    call function_2d_set ( 'SET', i )
!
!  Try an increasing number of subdivisions.
!
    do nsub = 1, 5

      call triangle_sub ( function_2d, xval, yval, nsub, order, xtab, ytab,  &
        weight, result )

      write ( *, '(2x,a7,i4, f14.8)' ) function_2d_name(i), nsub, result
 
    end do
  
  end do

  deallocate ( xtab )
  deallocate ( ytab )
  deallocate ( weight )

  return
end
subroutine test42

!*****************************************************************************80
!
!! TEST42 tests TRIANGLE_UNIT_SET and TRIANGLE_UNIT_SUM.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: rule_max = 20

  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) result(rule_max)
  integer   ( kind = 4 ) rule
  integer   ( kind = 4 ) triangle_unit_size
  real    ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real    ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real    ( kind = 8 ), allocatable, dimension ( : ) :: ytab

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST42'
  write ( *, '(a)' ) '  TRIANGLE_UNIT_SET sets up a quadrature '
  write ( *, '(a)' ) '    in the unit triangle,'
  write ( *, '(a)' ) '  TRIANGLE_UNIT_SUM applies it.'
  write ( *, '(a)' ) ' '

  do ilo = 1, rule_max, 5

    ihi = min ( ilo + 4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a,5(i7,7x))' ) 'Rule:    ', ( rule, rule = ilo, ihi )
    write ( *, '(2x,a)' ) 'Function '
    write ( *, '(a)' ) ' '

    num = function_2d_num ( )

    do i = 1, num

      call function_2d_set ( 'SET', i )

      do rule = ilo, ihi

        order = triangle_unit_size ( rule )

        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( weight(1:order) )

        call triangle_unit_set ( rule, order, xtab, ytab, weight )
 
        call triangle_unit_sum ( function_2d, order, xtab, ytab, weight, &
          result(rule) )

        deallocate ( xtab )
        deallocate ( ytab )
        deallocate ( weight )

      end do

      write ( *, '(2x,a7,5f14.8)' ) function_2d_name(i), result(ilo:ihi)

    end do

  end do

  return
end
subroutine test425

!*****************************************************************************80
!
!! TEST425 tests TRIANGLE_UNIT_SET.
!
!  Modified:
!
!    08 April 20008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) a
  integer ( kind = 4 ) b
  real    ( kind = 8 ) coef
  real    ( kind = 8 ) err
  real    ( kind = 8 ) exact
  integer ( kind = 4 ) i
  integer ( kind = 4 ) order
  real    ( kind = 8 ) quad
  integer ( kind = 4 ) rule
  integer ( kind = 4 ), parameter :: rule_max = 20
  integer ( kind = 4 ) triangle_unit_size
  real    ( kind = 8 ) value
  real    ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real    ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real    ( kind = 8 ), allocatable, dimension ( : ) :: ytab

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST425'
  write ( *, '(a)' ) '  TRIANGLE_UNIT_SET sets up a quadrature '
  write ( *, '(a)' ) '    in the unit triangle,'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Estimate integral of X^A * Y^B.'

  do a = 0, 10

    do b = 0, 10 - a

      coef = real ( a + b + 2, kind = 8 ) * real ( a + b + 1, kind = 8 )
      do i = 1, b
        coef = coef * real ( a + i, kind = 8 ) / real ( i, kind = 8 )
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a,i8,a,i8)' ) '  A = ', a, '  B = ', b
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Rule       QUAD           ERROR'
      write ( *, '(a)' ) ' '

      do rule = 1, rule_max

        order = triangle_unit_size ( rule )

        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( weight(1:order) )

        call triangle_unit_set ( rule, order, xtab, ytab, weight )
 
        quad = 0.0D+00

        do i = 1, order

          if ( a == 0 .and. b == 0 ) then
            value = coef
          else if ( a == 0 .and. b /= 0 ) then
            value = coef              * ytab(i)**b
          else if ( a /= 0 .and. b == 0 ) then
            value = coef * xtab(i)**a 
          else if ( a /= 0 .and. b /= 0 ) then
            value = coef * xtab(i)**a * ytab(i)**b
          end if

          quad = quad + 0.5D+00 * weight(i) * value

        end do

        exact = 1.0D+00
        err = abs ( exact - quad )

        write ( *, '(2x,i4,2x,g14.6,2x,f14.8)' ) rule, quad, err

        deallocate ( xtab )
        deallocate ( ytab )
        deallocate ( weight )

      end do

    end do

  end do

  return
end
subroutine test43

!*****************************************************************************80
!
!! TEST43 tests TRIANGLE_UNIT_PRODUCT_SET and TRIANGLE_UNIT_SUM.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: rule_max = 8

  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) result(rule_max)
  integer   ( kind = 4 ) rule
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST43'
  write ( *, '(a)' ) '  TRIANGLE_UNIT_PRODUCT_SET sets up a product quadrature'
  write ( *, '(a)' ) '    rule in the unit triangle,'
  write ( *, '(a)' ) '  TRIANGLE_UNIT_SUM applies it.'
  write ( *, '(a)' ) ' '

  do ilo = 1, rule_max, 5

    ihi = min ( ilo +  4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a,5i6)' ) 'Rule Order: ', ( rule, rule = ilo, ihi )
    write ( *, '(2x,a)' ) 'Function '
    write ( *, '(a)' ) ' '

    num = function_2d_num ( )

    do i = 1, num

      call function_2d_set ( 'SET', i )

      do rule = ilo, ihi

        call triangle_unit_product_size ( rule, order )

        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( weight(1:order) )

        call triangle_unit_product_set ( rule, order, xtab, ytab, weight )

        call triangle_unit_sum ( function_2d, order, xtab, ytab, weight, &
          result(rule) )

        deallocate ( xtab )
        deallocate ( ytab )
        deallocate ( weight )

      end do

      write ( *, '(2x,a7,5f14.8)' ) function_2d_name(i), result(ilo:ihi)

    end do

  end do

  return
end
subroutine test44

!*****************************************************************************80
!
!! TEST44 tests TRIANGLE_UNIT_SET and TRIANGLE_SUM.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ), parameter :: rule_max = 20

  real      ( kind = 8 ), external :: function_2d
  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) function_2d_num
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) ihi
  integer   ( kind = 4 ) ilo
  integer   ( kind = 4 ) num
  integer   ( kind = 4 ) order
  real      ( kind = 8 ) result(rule_max)
  integer   ( kind = 4 ) rule
  integer   ( kind = 4 ) triangle_unit_size
  real      ( kind = 8 ), allocatable, dimension ( : ) :: weight
  real      ( kind = 8 ), allocatable, dimension ( : ) :: xtab
  real      ( kind = 8 ), dimension ( 3 ) :: xval = (/ &
    1.0D+00, 3.0D+00, 1.0D+00 /)
  real      ( kind = 8 ), allocatable, dimension ( : ) :: ytab
  real      ( kind = 8 ), dimension ( 3 ) :: yval = (/ &
    1.0D+00, 1.0D+00, 4.0D+00 /)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST44'
  write ( *, '(a)' ) '  TRIANGLE_UNIT_SET sets up quadrature '
  write ( *, '(a)' ) '    in the unit triangle,'
  write ( *, '(a)' ) '  TRIANGLE_SUM applies it to an arbitrary triangle.'
  write ( *, '(a)' ) ' '

  do ilo = 1, rule_max, 5

    ihi = min ( ilo + 4, rule_max )

    write ( *, '(a)' ) ' '
    write ( *, '(2x,a,5i6)' ) 'Rule:    ', ( rule, rule = ilo, ihi )
    write ( *, '(2x,a)' ) 'Function '
    write ( *, '(a)' ) ' '

    num = function_2d_num ( )

    do i = 1, num

      call function_2d_set ( 'SET', i )

      do rule = ilo, ihi

        order = triangle_unit_size ( rule )

        allocate ( xtab(1:order) )
        allocate ( ytab(1:order) )
        allocate ( weight(1:order) )

        call triangle_unit_set ( rule, order, xtab, ytab, weight )
 
        call triangle_sum ( function_2d, xval, yval, order, xtab, ytab, weight, &
          result(rule) )

        deallocate ( xtab )
        deallocate ( ytab )
        deallocate ( weight )

      end do

      write ( *, '(2x,a7,5f14.8)' ) function_2d_name(i), result(ilo:ihi)

    end do

  end do

  return
end
subroutine setsim ( n, v )

!*****************************************************************************80
!
!! SETSIM defines a unit simplex.
!
!  Modified:
!
!    06 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  real    ( kind = 8 ) v(n,n+1)

  v(1:n,1:n+1) = 0.0D+00

  do i = 1, n
    v(i,i+1) = 1.0D+00
  end do
 
  return
end
subroutine test45

!*****************************************************************************80
!
!! TEST45 tests TORUS_1.
!
!  Modified:
!
!    05 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ) j2
  integer   ( kind = 4 ) n
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result(5)
  real      ( kind = 8 ) r1
  real      ( kind = 8 ) r2
  real      ( kind = 8 ) torus_area_3d

  r1 = 0.5D+00
  r2 = 1.0D+00
  n = 10
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST45'
  write ( *, '(a)' ) '  TORUS_1 approximates integrals on a torus.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The degree N will be varied.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Inner radius = ', r1
  write ( *, '(a,g14.6)' ) '  Outer radius = ', r2
  write ( *, '(a,g14.6)' ) '  Area = ', torus_area_3d ( r1, r2 )
  write ( *, '(a)' ) ' '
  write ( *, '(a,5i5)' ) '    F(X)  ', ( 2**j, j = 0, 8, 2 )
  write ( *, '(a)' ) ' '
 
  num = function_3d_num ( )

  do i = 1, num

    call function_3d_set ( 'SET', i )

    do j = 1, 5

      j2 = 2 * ( j - 1 )
      n = 2**j2
      call torus_1 ( function_3d, r1, r2, n, result(j) )

    end do

    write ( *, '(2x,a7,5f14.8)' ) function_3d_name(i), result(1:5)

  end do
 
  return
end
subroutine test46

!*****************************************************************************80
!
!! TEST46 tests TORUS_5S2, TORUS_6S2 and TORUS_14S.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) result3
  real      ( kind = 8 ) r1
  real      ( kind = 8 ) r2
  real      ( kind = 8 ) torus_volume_3d

  r1 = 0.5D+00
  r2 = 1.0D+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST46'
  write ( *, '(a)' ) '  For the interior of a torus,'
  write ( *, '(a)' ) '  TORUS_5S2,'
  write ( *, '(a)' ) '  TORUS_6S2, and'
  write ( *, '(a)' ) '  TORUS_5S2 approximate integrals.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Inner radius = ', r1
  write ( *, '(a,g14.6)' ) '  Outer radius = ', r2
  write ( *, '(a,g14.6)' ) '  Volume = ', torus_volume_3d ( r1, r2 )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    Rule:        #5S2          #6S2          #14S'
  write ( *, '(a)' ) '    F(X)'
  write ( *, '(a)' ) ' '
 
  num = function_3d_num ( )

  do i = 1, num

    call function_3d_set ( 'SET', i )

    call torus_5s2 ( function_3d, r1, r2, result1 )
    call torus_6s2 ( function_3d, r1, r2, result2 )
    call torus_14s ( function_3d, r1, r2, result3 )

    write ( *, '(2x,a7,3f14.8)' ) function_3d_name(i), result1, result2, result3

  end do
 
  return
end
subroutine test47

!*****************************************************************************80
!
!! TEST47 tests TORUS_SQUARE_5C2 and TORUS_SQUARE_14C.
!
!  Modified:
!
!    02 March 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  character ( len = 7 ) function_3d_name
  integer   ( kind = 4 ) function_3d_num
  real      ( kind = 8 ), external :: function_3d
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) num
  real      ( kind = 8 ) result1
  real      ( kind = 8 ) result2
  real      ( kind = 8 ) r1
  real      ( kind = 8 ) r2
  real      ( kind = 8 ) torus_square_volume_3d

  r1 = 1.0D+00
  r2 = 0.125D+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST47'
  write ( *, '(a)' ) '  For integrals inside a torus with square cross-section:'
  write ( *, '(a)' ) '  TORUS_SQUARE_5C2 approximates the integral;'
  write ( *, '(a)' ) '  TORUS_SQUARE_14C approximates the integral.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Inner radius = ', r1
  write ( *, '(a,g14.6)' ) '  Outer radius = ', r2
  write ( *, '(a,g14.6)' ) '  Volume = ', torus_square_volume_3d ( r1, r2 )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    F(X)    5C2           14C'
  write ( *, '(a)' ) ' '
 
  num = function_3d_num ( )

  do i = 1, num

    call function_3d_set ( 'SET', i )

    call torus_square_5c2 ( function_3d, r1, r2, result1 )
    call torus_square_14c ( function_3d, r1, r2, result2 )

    write ( *, '(2x,a7,2f14.8)' ) function_3d_name(i), result1, result2

  end do
 
  return
end
subroutine test48

!*****************************************************************************80
!
!! TEST48 tests TVEC_EVEN, TVEC_EVEN2 and TVEC_EVEN3.
!
!  Modified:
!
!    01 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) nt
  real    ( kind = 8 ), allocatable, dimension ( : ) :: t

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST48'
  write ( *, '(a)' ) '  For evenly spaced angles between 0 and 2*PI:'
  write ( *, '(a)' ) '  TVEC_EVEN'
  write ( *, '(a)' ) '  TVEC_EVEN2'
  write ( *, '(a)' ) '  TVEC_EVEN3'

  nt = 4
  allocate ( t(1:nt) )
  call tvec_even ( nt, t )
  call r8vec_print ( nt, t, '  TVEC_EVEN:' )
  deallocate ( t )

  nt = 4
  allocate ( t(1:nt) )
  call tvec_even2 ( nt, t )
  call r8vec_print ( nt, t, '  TVEC_EVEN2:' )
  deallocate ( t )

  nt = 4
  allocate ( t(1:nt) )
  call tvec_even3 ( nt, t )
  call r8vec_print ( nt, t, '  TVEC_EVEN3:' )
  deallocate ( t )

  return
end
subroutine test49

!*****************************************************************************80
!
!! TEST49 tests TVEC_EVEN_BRACKET, TVEC_EVEN_BRACKET2 and TVEC_EVEN_BRACKET3.
!
!  Modified:
!
!    05 April 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) nt
  real    ( kind = 8 ), allocatable, dimension ( : ) :: t
  real    ( kind = 8 ) theta1
  real    ( kind = 8 ) theta2

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST49'
  write ( *, '(a)' ) '  For evenly spaced angles between THETA1 and THETA2:'
  write ( *, '(a)' ) '  TVEC_EVEN_BRACKET'
  write ( *, '(a)' ) '  TVEC_EVEN_BRACKET2.'
  write ( *, '(a)' ) '  TVEC_EVEN_BRACKET3.'

  theta1 = 30.0D+00
  theta2 = 90.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '    THETA1 = ', theta1
  write ( *, '(a,g14.6)' ) '    THETA2 = ', theta2

  nt = 4
  allocate ( t(1:nt) )
  call tvec_even_bracket ( nt, theta1, theta2, t )
  call r8vec_print ( nt, t, '  TVEC_EVEN_BRACKET' )
  deallocate ( t )

  nt = 5
  allocate ( t(1:nt) )
  call tvec_even_bracket2 ( nt, theta1, theta2, t )
  call r8vec_print ( nt, t, '  TVEC_EVEN_BRACKET2' )
  deallocate ( t )

  nt = 3
  allocate ( t(1:nt) )
  call tvec_even_bracket3 ( nt, theta1, theta2, t )
  call r8vec_print ( nt, t, '  TVEC_EVEN_BRACKET3' )
  deallocate ( t )

  return
end
function function_1d ( x )

!*****************************************************************************80
!
!! FUNCTION_1D evaluates a function F(X) of one variable.
!
!  Discussion:
!
!    The actual form of the function can be determined by calling FUNCSET.
!
!  Modified:
!
!    21 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the value of the variable.
!
!    Output, real ( kind = 8 ) FUNCTION_1D, the value of the function.
!
  implicit none

  real    ( kind = 8 ) function_1d
  integer ( kind = 4 ) ifunc
  real    ( kind = 8 ) value
  real    ( kind = 8 ) x

  call function_1d_set ( 'GET', ifunc )

  if ( ifunc == 1 ) then
    value = 1.0D+00
  else if ( ifunc == 2 ) then
    value = x
  else if ( ifunc == 3 ) then
    value = x**2
  else if ( ifunc == 4 ) then
    value = x**3
  else if ( ifunc == 5 ) then
    value = x**4
  else if ( ifunc == 6 ) then
    value = x**5
  else if ( ifunc == 7 ) then
    value = x**6
  else if ( ifunc == 8 ) then
    value = abs ( x )
  else if ( ifunc == 9 ) then
    value = sin ( x )
  else if ( ifunc == 10 ) then
    value = exp ( x )
  else if ( ifunc == 11 ) then
    value = 1.0D+00 / ( 1.0D+00 + abs ( x ) )
  else if ( ifunc == 12 ) then
    value = sqrt ( abs ( x ) )
  else
    value = 0.0D+00
  end if

  function_1d = value

  return
end
function function_1d_name ( ifunc )

!*****************************************************************************80
!
!! FUNCTION_1D_NAME returns the name of the current 1D function.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IFUNC, the index of the function.
!
!    Output, character ( len = 7 ) FUNCTION_1D_NAME, the name of the 1D function.
!
  implicit none

  character ( len = 7 ) function_1d_name
  integer   ( kind = 4 ) ifunc

  if ( ifunc == 1 ) then
    function_1d_name = '      1'
  else if ( ifunc == 2 ) then
    function_1d_name = '      X'
  else if ( ifunc == 3 ) then
    function_1d_name = '    X^2'
  else if ( ifunc == 4 ) then
    function_1d_name = '    X^3'
  else if ( ifunc == 5 ) then
    function_1d_name = '    X^4'
  else if ( ifunc == 6 ) then
    function_1d_name = '    X^5'
  else if ( ifunc == 7 ) then
    function_1d_name = '    X^6'
  else if ( ifunc == 8 ) then
    function_1d_name = '      R'
  else if ( ifunc == 9 ) then
    function_1d_name = ' SIN(X)'
  else if ( ifunc == 10 ) then
    function_1d_name = ' EXP(X)'
  else if ( ifunc == 11 ) then
    function_1d_name = '1/(1+R)'
  else if ( ifunc == 12 ) then
    function_1d_name = 'SQRT(R)'
  else
    function_1d_name = '???????'
  end if

  return
end
function function_1d_num ( )

!*****************************************************************************80
!
!! FUNCTION_1D_NUM returns the number of 1D functions.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer ( kind = 4 ) FUNCTION_1D_NUM, the number of 1D functions.
!
  implicit none

  integer   ( kind = 4 ) function_1d_num

  function_1d_num = 12

  return
end
subroutine function_1d_set ( action, i )

!*****************************************************************************80
!
!! FUNCTION_1D_SET sets or reports the index of the current 1D function.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION:
!    'GET', please return the index of the current function.
!    'SET', please set the function according to the input index.
!
!    Input/output, integer I.
!    If ACTION = 'SET', then I is input, and is the index of the desired
!    function.
!    If ACTION = 'GET', then I is output, and is the index of the current
!    function.
!
  implicit none

  character ( len = * ) action
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ), save :: ival = 0

  if ( action == 'SET' ) then
    ival = i
  else if ( action == 'GET' ) then
    i = ival
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FUNCTION_1D_SET - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop
  end if

  return
end
function function_2d ( x, y )

!*****************************************************************************80
!
!! FUNCTION_2D evaluates a function F(X,Y) of two variables.
!
!  Modified:
!
!    19 April 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, Y, the value of the variables.
!
!    Output, real ( kind = 8 ) FUNCTION_2D, the value of the function.
!
  implicit none

  real    ( kind = 8 ) function_2d
  integer ( kind = 4 ) ifunc
  real    ( kind = 8 ) value
  real    ( kind = 8 ) x
  real    ( kind = 8 ) y

  call function_2d_set ( 'GET', ifunc )

  if ( ifunc == 1 ) then
    value = 1.0D+00
  else if ( ifunc == 2 ) then
    value = x
  else if ( ifunc == 3 ) then
    value = x * x
  else if ( ifunc == 4 ) then
    value = x * x * x
  else if ( ifunc == 5 ) then
    value = x * x * x * x
  else if ( ifunc == 6 ) then
    value = x**5
  else if ( ifunc == 7 ) then
    value = x**6
  else if ( ifunc == 8 ) then
    value = sqrt ( x * x + y * y )
  else if ( ifunc == 9 ) then
    value = sin ( x )
  else if ( ifunc == 10 ) then
    value = exp ( x )
  else if ( ifunc == 11 ) then
    value = 1.0D+00 / ( 1.0D+00 + sqrt ( x * x + y * y ) )
  else if ( ifunc == 12 ) then
    value = sqrt ( sqrt ( x * x + y * y ) )
  else
    value = 0.0D+00
  end if

  function_2d = value

  return
end
function function_2d_name ( ifunc )

!*****************************************************************************80
!
!! FUNCTION_2D_NAME returns the name of the current 2D function.
!
!  Modified:
!
!    07 April 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IFUNC, the index of the function.
!
!    Output, character ( len = 7 ) FUNCTION_2D_NAME, the name of the 2D function.
!
  implicit none

  character ( len = 7 ) function_2d_name
  integer   ( kind = 4 ) ifunc

  if ( ifunc == 1 ) then
    function_2d_name = '      1'
  else if ( ifunc == 2 ) then
    function_2d_name = '      X'
  else if ( ifunc == 3 ) then
    function_2d_name = '    X^2'
  else if ( ifunc == 4 ) then
    function_2d_name = '    X^3'
  else if ( ifunc == 5 ) then
    function_2d_name = '    X^4'
  else if ( ifunc == 6 ) then
    function_2d_name = '    X^5'
  else if ( ifunc == 7 ) then
    function_2d_name = '    X^6'
  else if ( ifunc == 8 ) then
    function_2d_name = '      R'
  else if ( ifunc == 9 ) then
    function_2d_name = ' SIN(X)'
  else if ( ifunc == 10 ) then
    function_2d_name = ' EXP(X)'
  else if ( ifunc == 11 ) then
    function_2d_name = '1/(1+R)'
  else if ( ifunc == 12 ) then
    function_2d_name = 'SQRT(R)'
  else
    function_2d_name = '???????'
  end if

  return
end
function function_2d_num ( )

!*****************************************************************************80
!
!! FUNCTION_2D_NUM returns the number of 2D functions.
!
!  Modified:
!
!    06 April 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer ( kind = 4 ) FUNCTION_2D_NUM, the number of 2D functions.
!
  implicit none

  integer   ( kind = 4 ) function_2d_num

  function_2d_num = 12

  return
end
subroutine function_2d_set ( action, i )

!*****************************************************************************80
!
!! FUNCTION_2D_SET sets or reports the index of the current 2D function.
!
!  Modified:
!
!    08 April 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION:
!    'GET', please return the index of the current function.
!    'SET', please set the function according to the input index.
!
!    Input/output, integer I.
!    If ACTION = 'SET', then I is input, and is the index of the desired
!    function.
!    If ACTION = 'GET', then I is output, and is the index of the current
!    function.
!
  implicit none

  character ( len = * ) action
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ), save :: ival = 0

  if ( action == 'SET' ) then
    ival = i
  else if ( action == 'GET' ) then
    i = ival
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FUNCTION_2D_SET - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop
  end if

  return
end
function function_3d ( x, y, z )

!*****************************************************************************80
!
!! FUNCTION_3D evaluates a the current 3D function.
!
!  Discussion:
!
!    The actual form of the function can be determined by calling FUNCTION_3D_SET.
!
!  Modified:
!
!    18 March 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, Y, Z, the value of the variables.
!
!    Output, real ( kind = 8 ) FUNCTION_3D, the value of the function.
!
  implicit none

  real    ( kind = 8 ) function_3d
  integer ( kind = 4 ) ifunc
  real    ( kind = 8 ) value
  real    ( kind = 8 ) x
  real    ( kind = 8 ) y
  real    ( kind = 8 ) z

  call function_3d_set ( 'GET', ifunc )

  if ( ifunc == 1 ) then
    value = 1.0D+00

  else if ( ifunc == 2 ) then
    value = x
  else if ( ifunc == 3 ) then
    value = y
  else if ( ifunc == 4 ) then
    value = z

  else if ( ifunc == 5 ) then
    value = x * x
  else if ( ifunc == 6 ) then
    value = x * y
  else if ( ifunc == 7 ) then
    value = x * z
  else if ( ifunc == 8 ) then
    value = y * y
  else if ( ifunc == 9 ) then
    value = y * z
  else if ( ifunc == 10 ) then
    value = z * z

  else if ( ifunc == 11 ) then
    value = x * x * x
  else if ( ifunc == 12 ) then
    value = x * y * z
  else if ( ifunc == 13 ) then
    value = z * z * z

  else if ( ifunc == 14 ) then
    value = x * x * x * x
  else if ( ifunc == 15 ) then
    value = x * x * z * z
  else if ( ifunc == 16 ) then
    value = z * z * z * z

  else if ( ifunc == 17 ) then
    value = x * x * x * x * x
  else if ( ifunc == 18 ) then
    value = x**6
  else if ( ifunc == 19 ) then
    value = sqrt ( x * x + y * y + z * z )
  else if ( ifunc == 20 ) then
    value = sin ( x )
  else if ( ifunc == 21 ) then
    value = exp ( x )
  else if ( ifunc == 22 ) then
    value = 1.0D+00 / sqrt ( 1.0D+00 + x * x + y * y + z * z )
  else if ( ifunc == 23 ) then
    value = sqrt ( sqrt ( x * x + y * y + z * z ) )
  else
    value = 0.0D+00
  end if

  function_3d = value

  return
end
function function_3d_name ( ifunc )

!*****************************************************************************80
!
!! FUNCTION_3D_NAME returns the name of the current 3D function.
!
!  Modified:
!
!    18 March 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, i