program dutch_prb
!
!*******************************************************************************
!
!! DUTCH_PRB tests the DUTCH routines.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DUTCH_PRB'
  write ( *, '(a)' ) '  Tests for DUTCH computational geometry routines.'
  write ( *, '(a)' ) ' '

  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test07
  call test08
  call test09
  call test10
  call test11
  call test12
  call test13
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DUTCH_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 tests POINTS_CONVEX_HULL_CUBIC_2D.
!
  implicit none
!
  integer, parameter :: nu = 7
!
  integer nv
  real ux(nu)
  real uy(nu)
  integer v(nu)
  real vx(nu)
  real vy(nu)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  POINTS_CONVEX_HULL_CUBIC_2D computes the convex hull of'
  write ( *, '(a)' ) '  a set of N 2D points using an algorithm that is cubic in N.'

  ux(1) = 0.0E+00
  uy(1) = 0.0E+00
  ux(2) = 1.0E+00
  uy(2) = 2.0E+00
  ux(3) = 2.0E+00
  uy(3) = 0.0E+00
  ux(4) = 1.0E+00
  uy(4) = 1.0E+00
  ux(5) = 0.0E+00
  uy(5) = 2.0E+00
  ux(6) = 1.0E+00
  uy(6) = 3.0E+00
  ux(7) = 2.0E+00
  uy(7) = 2.0E+00

  call rvec2_print ( nu, ux, uy, '  Coordinates of the points:' )

  call points_convex_hull_cubic_2d ( nu, ux, uy, nv, v )

  vx(1:nv) = ux(v(1:nv))
  vy(1:nv) = uy(v(1:nv))

  call rvec2_print ( nv, vx, vy, '  Coordinates of the convex hull:' )

  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 tests POINTS_CONVEX_HULL_NLOGN_2D.
!
  implicit none
!
  integer, parameter :: nu = 7
!
  integer nv
  real ux(nu)
  real uy(nu)
  integer v(nu)
  real vx(nu)
  real vy(nu)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  POINTS_CONVEX_HULL_NLOGN_2D computes the convex hull of'
  write ( *, '(a)' ) '  a set of N 2D points using an algorithm that is NlogN in N.'

  ux(1) = 0.0E+00
  uy(1) = 0.0E+00
  ux(2) = 1.0E+00
  uy(2) = 2.0E+00
  ux(3) = 2.0E+00
  uy(3) = 0.0E+00
  ux(4) = 1.0E+00
  uy(4) = 1.0E+00
  ux(5) = 0.0E+00
  uy(5) = 2.0E+00
  ux(6) = 1.0E+00
  uy(6) = 3.0E+00
  ux(7) = 2.0E+00
  uy(7) = 2.0E+00

  call rvec2_print ( nu, ux, uy, '  Coordinates of the points:' )

  call points_convex_hull_nlogn_2d ( nu, ux, uy, nv, v )

  vx(1:nv) = ux(v(1:nv))
  vy(1:nv) = uy(v(1:nv))

  call rvec2_print ( nv, vx, vy, '  Coordinates of the convex hull:' )

  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 tests POINTS_CONVEX_HULL_NLOGH_2D.
!
  implicit none
!
  integer, parameter :: nu = 7
!
  integer nv
  real ux(nu)
  real uy(nu)
  integer v(nu)
  real vx(nu)
  real vy(nu)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  POINTS_CONVEX_HULL_NLOGH_2D computes the convex hull of'
  write ( *, '(a)' ) '  a set of N 2D points using an algorithm that is NlogH in N.'
  write ( *, '(a)' ) '  (H is the number of points on the convex hull.)'

  ux(1) = 0.0E+00
  uy(1) = 0.0E+00
  ux(2) = 1.0E+00
  uy(2) = 2.0E+00
  ux(3) = 2.0E+00
  uy(3) = 0.0E+00
  ux(4) = 1.0E+00
  uy(4) = 1.0E+00
  ux(5) = 0.0E+00
  uy(5) = 2.0E+00
  ux(6) = 1.0E+00
  uy(6) = 3.0E+00
  ux(7) = 2.0E+00
  uy(7) = 2.0E+00

  call rvec2_print ( nu, ux, uy, '  Coordinates of the points:' )

  call points_convex_hull_nlogh_2d ( nu, ux, uy, nv, v )

  vx(1:nv) = ux(v(1:nv))
  vy(1:nv) = uy(v(1:nv))

  call rvec2_print ( nv, vx, vy, '  Coordinates of the convex hull:' )

  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 tests POLYCON_MINKOWSKI_SUM_LINEAR.
!
  implicit none
!
  integer, parameter :: nu = 4
  integer, parameter :: nv = 3
!
  integer nw
  real ux(nu)
  real uy(nu)
  real vx(nv)
  real vy(nv)
  real wx(nu+nv+2)
  real wy(nu+nv+2)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  POLYCON_MINKOWSKI_SUM_LINEAR computes the Minkowski sum of'
  write ( *, '(a)' ) '  two convex polygons using a linear algorithm.'

  ux(1) = 0.0E+00
  uy(1) = 0.0E+00
  ux(2) = 2.0E+00
  uy(2) = 2.0E+00
  ux(3) = -1.0E+00
  uy(3) = 3.0E+00
  ux(4) = -2.0E+00
  uy(4) = 2.0E+00

  call rvec2_print ( nu, ux, uy, '  Coordinates of polygon U:' )

  vx(1) = 8.0E+00
  vy(1) = 2.0E+00
  vx(2) = 9.0E+00
  vy(2) = 5.0E+00
  vx(3) = 7.0E+00
  vy(3) = 4.0E+00

  call rvec2_print ( nv, vx, vy, '  Coordinates of polygon V:' )
  
  call polycon_minkowski_sum_linear ( nu, ux, uy, nv, vx, vy, nw, wx, wy )

  call rvec2_print ( nw, wx, wy, '  Coordinates of Minkowski sum polygon W:' )

  return
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 tests POLYCON_MINKOWSKI_SUM_LINEAR.
!
  implicit none
!
  integer, parameter :: nu = 4
  integer, parameter :: nv = 3
!
  integer nw
  real ux(nu)
  real uy(nu)
  real vx(nv)
  real vy(nv)
  real wx(nu+nv+2)
  real wy(nu+nv+2)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  POLYCON_MINKOWSKI_SUM_N2LOGN2 computes the Minkowski'
  write ( *, '(a)' ) '  sum of two convex polygons using a linear algorithm.'

  ux(1) = 0.0E+00
  uy(1) = 0.0E+00
  ux(2) = 2.0E+00
  uy(2) = 2.0E+00
  ux(3) = -1.0E+00
  uy(3) = 3.0E+00
  ux(4) = -2.0E+00
  uy(4) = 2.0E+00

  call rvec2_print ( nu, ux, uy, '  Coordinates of polygon U:' )

  vx(1) = 8.0E+00
  vy(1) = 2.0E+00
  vx(2) = 9.0E+00
  vy(2) = 5.0E+00
  vx(3) = 7.0E+00
  vy(3) = 4.0E+00

  call rvec2_print ( nv, vx, vy, '  Coordinates of polygon V:' )
  
  call polycon_minkowski_sum_n2logn2 ( nu, ux, uy, nv, vx, vy, nw, wx, wy )

  call rvec2_print ( nw, wx, wy, '  Coordinates of Minkowski sum polygon W:' )

  return
end
subroutine test06
!
!*******************************************************************************
!
!! TEST06 tests PERM_RANDOM.
!
  implicit none
!
  integer, parameter :: n = 4
!
  integer i
  integer p(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  PERM_RANDOM produces a random permutation;'
  write ( *, '(a,i6)' ) '  For this test, N = ', n
  write ( *, '(a)' ) ' '
 
  do i = 1, 5
    call perm_random ( n, p )
    call perm_print ( n, p, ' ' )
  end do
  
  return
end
subroutine test07
!
!*******************************************************************************
!
!! TEST07 tests POINTS_MINIDISC.
!
  implicit none
!
  integer, parameter :: n = 7
!
  real cx
  real cy
  real px(n)
  real py(n)
  real r
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  POINTS_MINIDISC computes the smallest circle that'
  write ( *, '(a)' ) '  contains a set of N 2D points.'

  px(1) = 0.0E+00
  py(1) = 0.0E+00
  px(2) = 1.0E+00
  py(2) = 2.0E+00
  px(3) = 2.0E+00
  py(3) = 0.0E+00
  px(4) = 1.0E+00
  py(4) = 1.0E+00
  px(5) = 0.0E+00
  py(5) = 2.0E+00
  px(6) = 1.0E+00
  py(6) = 3.0E+00
  px(7) = 2.0E+00
  py(7) = 2.0E+00

  call rvec2_print ( n, px, py, '  Coordinates of the points:' )

  call points_minidisc_2d ( n, px, py, r, cx, cy )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The enclosing circle has:'
  write ( *, '(a,g14.6)' ) '  Radius R =      ', r
  write ( *, '(a,2g14.6)' ) '  Center (CX,CY)= ', cx, cy

  return
end
subroutine test08
!
!*******************************************************************************
!
!! TEST08 tests LINE_SEG_VEC_INT_2D.
!
!  Expecting:
!
!    2, 4, 2.0, 1.0
!    2, 8, 2.0, 1.0
!    4, 6, 3.0, 2.0
!    4, 8, 2.0, 1.0
!    6, 8, 2.0, 3.0
!
  implicit none
!
  integer, parameter :: n = 8
!
  integer flag
  integer i
  integer j
  real xint
  real, dimension ( n ) :: x1 = &
    (/ 4.0E+00, 0.0E+00, 4.0E+00, 1.0E+00, 0.0E+00, 1.0E+00, 3.0E+00, 2.0E+00 /)
  real, dimension ( n ) :: x2 = &
    (/ 5.0E+00, 3.0E+00, 5.0E+00, 5.0E+00, 1.0E+00, 4.0E+00, 3.0E+00, 2.0E+00 /)
  real yint
  real, dimension ( n ) :: y1 = &
    (/ 0.0E+00, 1.0E+00, 2.0E+00, 0.0E+00, 3.0E+00, 4.0E+00, 4.0E+00, 4.0E+00 /)
  real, dimension ( n ) :: y2 = &
    (/ 0.0E+00, 1.0E+00, 3.0E+00, 4.0E+00, 2.0E+00, 1.0E+00, 3.0E+00, 0.0E+00 /)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  LINE_SEG_VEC_INT_2D finds the intersections of a set'
  write ( *, '(a)' ) '  of line segments.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I, X1(I), Y1(I), X2(I), Y2(I)'
  write ( *, '(a)' ) ' ' 
  do i = 1, n
    write ( *, '(i2,4f7.1)' ) i, x1(i), y1(i), x2(i), y2(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Seg#1, Seg#2, Intersection'
  write ( *, '(a)' ) ' '

  i = 0
  j = 0

  do

    call line_seg_vec_int_2d ( n, x1, y1, x2, y2, i, j, flag, xint, yint )

    if ( flag == 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  No more intersections.'
      exit
    end if

    write ( *, '(2i6,2g14.6)' ) i, j, xint, yint

  end do

  return
end
subroutine test09
!
!*******************************************************************************
!
!! TEST09 tests RECT_INT_2D.
!
  implicit none
!
  integer, parameter :: ntest = 8
!
  integer flag
  integer i
  real x1
  real x2
  real, dimension ( ntest ) :: x3 = &
    (/ 1.0E+00, 3.0E+00, -2.0E+00, 3.0E+00, 2.0E+00, 5.0E+00, 5.0E+00, 6.0E+00 /)
  real, dimension ( ntest ) :: x4 = &
    (/ 2.0E+00, 5.0E+00,  1.0E+00, 9.0E+00, 8.0E+00, 8.0E+00, 8.0E+00, 8.0E+00 /)
  real x5
  real x6
  real y1
  real y2
  real, dimension ( ntest ) :: y3 = (/ &
    2.0E+00, 2.0E+00, 2.0E+00, -2.0E+00, -3.0E+00, &
    1.0E+00, 5.0E+00, 7.0E+00 /)
  real, dimension ( ntest ) :: y4 = (/ &
    4.0E+00, 5.0E+00, 4.0E+00,  1.0E+00,  8.0E+00, &
    3.0E+00, 8.0E+00, 9.0E+00 /)
  real y5
  real y6
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  RECT_INT_2D finds the intersections of two rectangles.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  For all tests, rectangle #1 is (0,0), (5,5).'
  write ( *, '(a)' ) ' '

  x1 = 0.0E+00
  y1 = 0.0E+00

  x2 = 5.0E+00
  y2 = 5.0E+00

  do i = 1, ntest

    call rect_int_2d ( x1, y1, x2, y2, x3(i), y3(i), x4(i), y4(i), &
      flag, x5, y5, x6, y6 )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Second rectangle:'
    write ( *, '(a)' ) ' '
    write ( *, '(2f8.4)' ) x3(i), y3(i)
    write ( *, '(2f8.4)' ) x4(i), y4(i)
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Intersection rectangle:'
    write ( *, '(a)' ) ' '

    if ( flag == 0 ) then
      write ( *, '(a)' ) '  EMPTY'
    else
      write ( *, '(2f8.4)' ) x5, y5
      write ( *, '(2f8.4)' ) x6, y6
    end if

  end do

  return
end
subroutine test10
!
!*******************************************************************************
!
!! TEST10 tests POLY_REORDER_NODES.
!
  implicit none
!
  integer, parameter :: npoly = 7
  integer, parameter :: nxy = 8
!
  integer, dimension ( npoly ) :: poly = (/ 7, 8, 1, 6, 2, 5, 4 /)
  real, dimension ( nxy ) :: x = &
    (/ 4.0E+00, 0.0E+00, 0.0E+00, 1.0E+00, 2.0E+00, 2.0E+00, 3.0E+00, 5.0E+00 /)
  real, dimension ( nxy ) :: y = &
    (/ 3.0E+00, 3.0E+00, 0.0E+00, 0.0E+00, 2.0E+00, 5.0E+00, 1.0E+00, 0.0E+00 /)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  POLY_REORDER_NODES reorders the nodes of a polygon'
  write ( *, '(a)' ) '  so that node 1 is leftmost (and lowest, in ties).'

  call rvec2_print ( nxy, x, y, '  The nodes:' )

  call ivec_print ( npoly, poly, '  The sequence of nodes:' )

  call poly_reorder_nodes ( nxy, x, y, npoly, poly )

  call ivec_print ( npoly, poly, '  The reordered sequence of nodes:' )

  return
end
subroutine test11
!
!*******************************************************************************
!
!! TEST11 tests POLY_TRIANGULATE_2D.
!
  implicit none
!
  integer, parameter :: n = 13
!
  integer i
  integer triang(3,n-2)
  real, dimension ( n ) :: x = (/ &
    7.0E+00, 7.0E+00, 9.0E+00,  4.0E+00, 4.0E+00, &
    6.0E+00, 6.0E+00, 4.0E+00,  2.0E+00, 2.0E+00, &
    0.0E+00, 2.0E+00, 4.0E+00 /)
  real, dimension ( n ) :: y = (/ &   
    0.0E+00, 4.0E+00, 5.0E+00, 10.0E+00, 4.0E+00, &
    6.0E+00, 2.0E+00, 2.0E+00,  4.0E+00, 2.0E+00, &
    2.0E+00, 0.0E+00, 1.0E+00 /)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  POLY_TRIANGULATE_2D triangulates a polygon.'

  call rvec2_print ( n, x, y, '  The nodes of the polygon:' )

  call poly_triangulate_2d ( n, x, y, triang )

  call imat_print ( 3, 3, n-2, triang, '  The triangulation:' )

  return
end
subroutine test12
!
!*******************************************************************************
!
!! TEST12 tests TRIANGULATE_TRICOLOR.
!
  implicit none
!
  integer, parameter :: n = 13
!
  integer color(n)
  integer i
  integer, dimension ( 3, n-2 ) :: triang =  reshape ( (/ &
    11,    12,    10, &
    12,    13,    10, &
    10,    13,     9, &
     9,    13,     8, &
    13,     1,     8, &
     8,     1,     7, &
     5,     6,     4, &
     4,     6,     3, &
     7,     1,     6, &
     6,     2,     3, &
     6,     1,     2 /), (/ 3, n-2 /) )
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  TRIANGULATE_TRICOLOR tricolors a triangulation.'

  call imat_print ( 3, 3, n-2, triang, '  The triangulation:' )

  call triangulate_tricolor ( n, triang, color )

  call ivec_print ( n, color, '  The node coloring' )

  return
end
subroutine test13
!
!*******************************************************************************
!
!! TEST13 tests TRIANGULATION_BOUNDARY_COUNT.
!
  implicit none
!
  integer bound_num
  integer, parameter :: point_num = 13
  integer, parameter :: tri_num = 16
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST13'
  write ( *, '(a)' ) '  TRIANGULATION_BOUNDARY_COUNT determines the number of'
  write ( *, '(a)' ) '  edges that lie on the convex hull of a region that has'
  write ( *, '(a)' ) '  been triangulated.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of points =         ', point_num
  write ( *, '(a,i6)' ) '  Number of triangles =      ', tri_num

  call triangulation_boundary_count ( point_num, tri_num, bound_num )

  write ( *, '(a,i6)' ) '  Number of boundary edges = ', bound_num

  return
end