program main

!*****************************************************************************80
!
!! test_int_2d_test() tests test_int_2d().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 January 2009
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'test_int_2d_test():'
  write ( *, '(a)' ) '  Fortran90 version'
  write ( *, '(a)' ) '  Test TEST_INT_2D().'

  call test01 ( )
  call test02 ( )
  call test03 ( )
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST_INT_2D_TEST():'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop 0
end
subroutine test01 ( )

!*****************************************************************************80
!
!! TEST01 applies a Monte Carlo rule.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 January 2009
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) a(2)
  real ( kind = rk ) b(2)
  integer dim
  real ( kind = rk ) error
  real ( kind = rk ) exact
  real ( kind = rk ), allocatable :: fx(:)
  integer i
  integer n
  integer problem
  integer problem_num
  real ( kind = rk ) quad
  integer seed
  real ( kind = rk ) volume
  real ( kind = rk ), allocatable :: x(:,:)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  Use a Monte Carlo rule.'
  write ( *, '(a)' ) '  Repeatedly multiply the number of points by 4.'

  call p00_problem_num ( problem_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) &
    '   Problem      Points       Approx         Error'

  do problem = 1, problem_num

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

    n = 1

    do i = 1, 12

      seed = 123456789

      allocate ( x(2,n) )
      allocate ( fx(n) )

      call random_number ( harvest = x(1:2,1:n) )

      call p00_lim ( problem, a, b )

      do dim = 1, 2
        x(dim,1:n) = ( 1.0D+00 - x(dim,1:n) ) * a(dim) &
                   +             x(dim,1:n)   * b(dim)
      end do

      volume = product ( b(1:2) - a(1:2) )

      call p00_fun ( problem, n, x, fx )

      quad = volume * sum ( fx(1:n) ) / real ( n, kind = rk )

      call p00_exact ( problem, exact )

      error = abs ( quad - exact )

      write ( *, '(2x,i8,2x,i10,2x,g14.6,2x,g14.6)' ) &
        problem, n, quad, error

      deallocate ( fx )
      deallocate ( x )

      n = n * 4

    end do

    write ( *, '(2x,i8,2x,a10,2x,g14.6)' ) &
      problem, '     Exact', exact

  end do

  return
end
subroutine test02 ( )

!*****************************************************************************80
!
!! TEST02 applies a product of composite midpoint rules.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 July 2010
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) a(2)
  real ( kind = rk ) b(2)
  real ( kind = rk ) error
  real ( kind = rk ) exact
  real ( kind = rk ), allocatable :: fx(:)
  integer i
  integer ix
  integer iy
  integer k
  integer n
  integer nx
  integer ny
  integer problem
  integer problem_num
  real ( kind = rk ) quad
  real ( kind = rk ) volume
  real ( kind = rk ), allocatable :: x(:,:)
  real ( kind = rk ) xval
  real ( kind = rk ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Use a product of composite midpoint rules.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Repeatedly multiply the number of points by 4.'

  call p00_problem_num ( problem_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) &
    '   Problem      Points       Approx         Error'

  do problem = 1, problem_num

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

    nx = 1
    ny = 1

    do i = 1, 12

      n = nx * ny

      allocate ( x(2,n) )
      allocate ( fx(n) )

      call p00_lim ( problem, a, b )

      k = 0

      do ix = 1, nx

        xval = ( real ( 2 * nx - 2 * ix + 1, kind = rk ) * a(1)   &
               + real (          2 * ix - 1, kind = rk ) * b(1) ) &
               / real ( 2 * nx,              kind = rk )

        do iy = 1, ny

          yval = ( real ( 2 * ny - 2 * iy + 1, kind = rk ) * a(2)   &
                 + real (          2 * iy - 1, kind = rk ) * b(2) ) &
                 / real ( 2 * ny,              kind = rk )

          k = k + 1
          x(1,k) = xval
          x(2,k) = yval

        end do

      end do

      volume = product ( b(1:2) - a(1:2) )

      call p00_fun ( problem, n, x, fx )

      quad = volume * sum ( fx(1:n) ) / real ( n, kind = rk )

      call p00_exact ( problem, exact )

      error = abs ( quad - exact )

      write ( *, '(2x,i8,2x,i10,2x,g14.6,2x,g14.6)' ) &
        problem, n, quad, error

      deallocate ( fx )
      deallocate ( x )

      nx = nx * 2
      ny = ny * 2

    end do

    write ( *, '(2x,i8,2x,a10,2x,g14.6)' ) &
      problem, '     Exact', exact

  end do

  return
end
subroutine test03 ( )

!*****************************************************************************80
!
!! TEST03 applies a product of Gauss-Legendre rules.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 September 2011
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) a(2)
  real ( kind = rk ) b(2)
  real ( kind = rk ) error
  real ( kind = rk ) exact
  real ( kind = rk ), allocatable :: fxy(:)
  integer i
  integer ix
  integer iy
  integer k
  integer nx
  integer nxy
  integer ny
  integer problem
  integer problem_num
  real ( kind = rk ) quad
  real ( kind = rk ) volume
  real ( kind = rk ), allocatable :: w(:)
  real ( kind = rk ), allocatable :: wxy(:)
  real ( kind = rk ), allocatable :: x(:)
  real ( kind = rk ), allocatable :: xy(:,:)
  real ( kind = rk ) xval
  real ( kind = rk ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  Use a product of Gauss-Legendre rules.'
  write ( *, '(a)' ) '  The 1D rules essentially double in order.'

  call p00_problem_num ( problem_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) &
    '   Problem      Points       Approx         Error'

  do problem = 1, problem_num

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

    nx = 1
    ny = 1

    do i = 1, 8

      allocate ( x(1:nx) )
      allocate ( w(1:nx) )

      call legendre_set ( nx, x, w )

      nxy = nx * ny

      allocate ( wxy(nxy) )
      allocate ( xy(2,nxy) )
      allocate ( fxy(nxy) )

      call p00_lim ( problem, a, b )

      k = 0

      do ix = 1, nx

        xval = ( ( 1.0D+00 + x(ix) ) * a(1)   &
               + ( 1.0D+00 - x(ix) ) * b(1) ) &
               /   2.0D+00

        do iy = 1, ny

          yval = ( ( 1.0D+00 + x(iy) ) * a(2)   &
                 + ( 1.0D+00 - x(iy) ) * b(2) ) &
                 /   2.0D+00

          k = k + 1
          xy(1,k) = xval
          xy(2,k) = yval
          wxy(k) = w(ix) * w(iy)

        end do

      end do

      volume = product ( b(1:2) - a(1:2) )

      call p00_fun ( problem, nxy, xy, fxy )

      quad = volume * dot_product ( wxy(1:nxy), fxy(1:nxy) ) / 4.0D+00

      call p00_exact ( problem, exact )

      error = abs ( quad - exact )

      write ( *, '(2x,i8,2x,i10,2x,g14.6,2x,g14.6)' ) &
        problem, nxy, quad, error

      deallocate ( fxy )
      deallocate ( w )
      deallocate ( wxy )
      deallocate ( x )
      deallocate ( xy )

      nx = 2 * nx + 1
      ny = nx

    end do

    write ( *, '(2x,i8,2x,a10,2x,g14.6)' ) &
      problem, '     Exact', exact

  end do

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! timestamp() prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end