function angle ( xa, ya, xb, yb, xc, yc )

!*****************************************************************************80
!
!! ANGLE computes the interior angle at a vertex defined by 3 points.
!
!  Discussion:
!
!    ANGLE computes the interior angle, in radians, at vertex
!    (XB,YB) of the chain formed by the directed edges from
!    (XA,YA) to (XB,YB) to (XC,YC).  The interior is to the
!    left of the two directed edges.
!
!  Modified:
!
!    17 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) XA, YA, XB, YB, XC, YC, the coordinates of the 
!    vertices.
!
!    Output, real ( kind = 8 ) ANGLE, the interior angle formed by
!    the vertex, in radians, between 0 and 2*PI.
!
  implicit none

  real ( kind = 8 ) angle
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) t
  real ( kind = 8 ) x1
  real ( kind = 8 ) x2
  real ( kind = 8 ) xa
  real ( kind = 8 ) xb
  real ( kind = 8 ) xc
  real ( kind = 8 ) y1
  real ( kind = 8 ) y2
  real ( kind = 8 ) ya
  real ( kind = 8 ) yb
  real ( kind = 8 ) yc

  x1 = xa - xb
  y1 = ya - yb
  x2 = xc - xb
  y2 = yc - yb

  t = sqrt ( ( x1 * x1 + y1 * y1 ) * ( x2 * x2 + y2 * y2 ) )

  if ( t == 0.0D+00 ) then
    angle = pi
    return
  end if

  t = ( x1 * x2 + y1 * y2 ) / t

  if ( t < -1.0D+00 ) then
    t = -1.0D+00
  else if ( 1.0D+00 < t ) then
    t = 1.0D+00
  end if

  angle = acos ( t )

  if ( x2 * y1 - y2 * x1 < 0.0D+00 ) then
    angle = 2.0D+00 * pi - angle
  end if

  return
end
function areapg ( nvrt, xc, yc )

!*****************************************************************************80
!
!! AREAPG computes twice the signed area of a simple polygon.
!
!  Modified:
!
!    13 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer NVRT, the number of vertices on the boundary of 
!    the polygon.  N must be at least 3.
!
!    Input, real ( kind = 8 ) XC(NVRT), YC(NVRT), the X and Y coordinates 
!    of the vertices.
!
!    Output, real ( kind = 8 ) AREAPG, twice the signed area of the polygon,
!    which will be positive if the vertices were listed in counter clockwise
!    order, and negative otherwise.
!
  implicit none

  integer nvrt

  real ( kind = 8 ) areapg
  integer i
  real ( kind = 8 ) sum2
  real ( kind = 8 ) xc(nvrt)
  real ( kind = 8 ) yc(nvrt)

  sum2 = xc(1) * ( yc(2) - yc(nvrt) )

  do i = 2, nvrt-1
    sum2 = sum2 + xc(i) * ( yc(i+1) - yc(i-1) )
  end do

  sum2 = sum2 + xc(nvrt) * ( yc(1) - yc(nvrt-1) )

  areapg = sum2

  return
end
function areatr ( xa, ya, xb, yb, xc, yc )

!*****************************************************************************80
!
!! AREATR computes twice the signed area of a triangle.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) XA, YA, XB, YB, XC, YC, the coordinates of the 
!    vertices.
!
!    Output, real ( kind = 8 ) AREATR, twice the signed area of the triangle.
!    This will be positive if the vertices are listed in counter clockwise 
!    order.
!
  implicit none

  real ( kind = 8 ) areatr 
  real ( kind = 8 ) xa
  real ( kind = 8 ) xb
  real ( kind = 8 ) xc
  real ( kind = 8 ) ya
  real ( kind = 8 ) yb
  real ( kind = 8 ) yc

  areatr = ( xb - xa ) * ( yc - ya ) - ( xc - xa ) * ( yb - ya )

  return
end
subroutine bedgmv ( nvc, npolg, nvert, maxvc, h, vcl, hvl, pvl, vstart, vnum, &
  ierror )

!*****************************************************************************80
!
!! BEDGMV generates boundary edge mesh vertices.
!
!  Purpose: 
!
!    Generate mesh vertices on boundary of convex polygons
!    of decomposition with spacing determined by H array.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input/output, integer NVC, the number of coordinates or positions used 
!    in VCL array.
!
!    Input, integer NPOLG, the number of polygons or positions used 
!    in HVL array.
!
!    Input, integer NVERT, the number of vertices or positions used in 
!    PVL array.
!
!    Input, integer MAXVC, the maximum size available for VCL array.
!
!    Input, real ( kind = 8 ) H(1:NPOLG), the spacing of mesh vertices for 
!    convex polygons.
!
!    Input/output, real ( kind = 8 ) VCL(1:2,1:NVC), the vertex coordinate list.
!
!    Input, integer HVL(1:NPOLG, the head vertex list.
!
!    Input, integer PVL(1:4,1:NVERT), the polygon vertex list.
!
!    Output, integer VSTART(1:NVERT), the start location in VCL for mesh 
!    vertices on each edge in PVL if there are any, else 0.
!
!    Output, integer VNUM(1:NVERT), the number of mesh vertices on interior 
!    of each edge in PVL; entry is negated if mesh vertices are listed in 
!    backward order in VCL.
!
!    Output, integer IERROR, is set to 3 on error.
!
  implicit none

  integer maxvc
  integer npolg
  integer nvert

  real ( kind = 8 ) dx
  real ( kind = 8 ) dy
  integer, parameter :: edgv = 4
  real ( kind = 8 ) h(npolg)
  real ( kind = 8 ) hh
  integer hvl(npolg)
  integer i
  integer ia
  integer ierror
  integer j
  integer k
  integer l
  real ( kind = 8 ) leng
  integer, parameter :: loc = 1
  integer m
  integer nvc
  integer, parameter :: polg = 2
  integer pvl(4,nvert)
  integer, parameter :: succ = 3
  integer u
  integer v
  integer vstart(nvert)
  integer vnum(nvert)
  real ( kind = 8 ) vcl(2,maxvc)
  real ( kind = 8 ) x
  real ( kind = 8 ) y

  ierror = 0
  vstart(1:nvert) = -1

  do k = 1, npolg

    i = hvl(k)

    do

      j = pvl(succ,i)

      if ( vstart(i) == -1 ) then

        u = pvl(loc,i)
        v = pvl(loc,j)
        x = vcl(1,u)
        y = vcl(2,u)
        leng = sqrt ( ( vcl(1,v) - x )**2 + ( vcl(2,v) - y )**2 )
        ia = pvl(edgv,i)

        if ( ia <= 0 ) then
          hh = h(k)
        else
          hh = sqrt ( h(k) * h(pvl(polg,ia)) )
        end if

        if ( hh == 0.0D+00 ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'BEDGMV - Fatal error!'
          write ( *, '(a)' ) '  HH = 0.'
          stop
        end if

        l = int ( leng / hh )

        if ( real ( l, kind = 8 ) / real ( 2 * l + 1, kind = 8 ) &
          < ( leng / hh ) - real ( l, kind = 8 ) ) then
          l = l + 1
        end if

        if ( l <= 1 ) then

          vstart(i) = 0
          vnum(i) = 0

        else

          dx = ( vcl(1,v) - x ) / real ( l, kind = 8 )
          dy = ( vcl(2,v) - y ) / real ( l, kind = 8 )
          l = l - 1

          if ( maxvc < nvc + l ) then
            write ( *, '(a)' ) ' '
            write ( *, '(a)' ) 'BEDGMV - Fatal error!'
            write ( *, '(a)' ) '  MAXVC < NVC + L.'
            ierror = 3
            return
          end if

          vstart(i) = nvc + 1
          vnum(i) = l

          do m = 1, l
            x = x + dx
            y = y + dy
            nvc = nvc + 1
            vcl(1,nvc) = x
            vcl(2,nvc) = y
          end do

        end if

        if ( 0 < ia ) then
          vstart(ia) = vstart(i)
          vnum(ia) = -vnum(i)
        end if

      end if

      i = j

      if ( i == hvl(k) ) then
        exit
      end if

    end do

  end do

  return
end
subroutine bnsrt2 ( binexp, n, a, map, bin, iwk )

!*****************************************************************************80
!
!! BNSRT2 bin sorts N points in 2D into increasing bin order.
!
!  Purpose: 
!
!    Use a bin sort to obtain the permutation of N 2D
!    double precision points so that points are in increasing bin
!    order, where the N points are assigned to about N**BINEXP bins.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer BINEXP, the exponent for the number of bins.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) A(2,*), the points to be binned.
!
!    Input/output, integer MAP(N); on input, the points of A with indices 
!    MAP(1), MAP(2), ..., MAP(N) are to be sorted.  On output, MAP has
!    been permuted so bin of MAP(1) <= bin of MAP(2) <= ... <= bin of MAP(N).
!
!    Workspace, integer BIN(N), used for bin numbers and permutation of 1 to N.
!
!    Workspace, integer IWK(N), used for copy of MAP array.
!
  implicit none

  integer n

  real ( kind = 8 ) a(2,*)
  integer bin(n)
  real ( kind = 8 ) binexp
  real ( kind = 8 ) dx
  real ( kind = 8 ) dy
  integer i
  integer iwk(n)
  integer j
  integer k
  integer l
  integer map(n)
  integer nside
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin

  nside = int ( real ( n, kind = 8 )**( binexp / 2.0D+00 ) + 0.5D+00 )

  if ( nside <= 1 ) then
    return
  end if

  xmin = a(1,map(1))
  ymin = a(2,map(1))
  xmax = xmin
  ymax = ymin
  do i = 2, n
    j = map(i)
    xmin = min ( xmin, a(1,j) )
    xmax = max ( xmax, a(1,j) )
    ymin = min ( ymin, a(2,j) )
    ymax = max ( ymax, a(2,j) )
  end do

  dx = 1.0001D+00 * ( xmax - xmin ) / real ( nside, kind = 8 )
  dy = 1.0001D+00 * ( ymax - ymin ) / real ( nside, kind = 8 )

  if ( dx == 0.0D+00 ) then
    dx = 1.0D+00
  end if

  if ( dy == 0.0D+00 ) then
    dy = 1.0D+00
  end if

  do i = 1, n
    j = map(i)
    iwk(i) = j
    map(i) = i
    k = int ( ( a(1,j) - xmin ) / dx )
    l = int ( ( a(2,j) - ymin ) / dy )
    if ( mod ( k, 2 ) == 0 ) then
      bin(i) = k * nside + l
    else
      bin(i) = ( k + 1 ) * nside - l - 1
    end if
  end do

  call ihpsrt ( 1, n, 1, bin, map )

  bin(1:n) = map(1:n)

  do i = 1, n
    map(i) = iwk(bin(i))
  end do

  return
end
function cmcirc ( x0, y0, x1, y1, x2, y2, x3, y3 )

!*****************************************************************************80
!
!! CMCIRC determines whether a point lies within a circle through 3 points.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X0, Y0, the coordinates of the point to
!    be tested.
!
!    Input, real ( kind = 8 ) X1, Y1, X2, Y2, X3, Y3, the coordinates of
!    three points that define a circle.
!
!    Output, integer CMCIRC, reports the test results:
!     2, if the three vertices are collinear,
!     1, if (X0,Y0) is inside the circle,
!     0, if (X0,Y0) is on the circle,
!    -1, if (X0,Y0) is outside the circle.
!
  real ( kind = 8 ) a11
  real ( kind = 8 ) a12
  real ( kind = 8 ) a21
  real ( kind = 8 ) a22
  real ( kind = 8 ) b1
  real ( kind = 8 ) b2
  integer cmcirc
  real ( kind = 8 ) det
  real ( kind = 8 ) diff
  real ( kind = 8 ) rsq
  real ( kind = 8 ) tol
  real ( kind = 8 ) tolabs
  real ( kind = 8 ) xc
  real ( kind = 8 ) yc
  real ( kind = 8 ) x0
  real ( kind = 8 ) x1
  real ( kind = 8 ) x2
  real ( kind = 8 ) x3
  real ( kind = 8 ) y0
  real ( kind = 8 ) y1
  real ( kind = 8 ) y2
  real ( kind = 8 ) y3

  tol = 100.0D+00 * epsilon ( tol )
  cmcirc = 2
  a11 = x2 - x1
  a12 = y2 - y1
  a21 = x3 - x1
  a22 = y3 - y1
  tolabs = tol * max ( abs ( a11), abs ( a12), abs ( a21), abs ( a22) )
  det = a11 * a22 - a21 * a12

  if ( abs ( det ) <= tolabs ) then
    return
  end if

  b1 = a11 * a11 + a12 * a12
  b2 = a21 * a21 + a22 * a22
  det = 2.0D+00 * det
  xc = ( b1 * a22 - b2 * a12 ) / det
  yc = ( b2 * a11 - b1 * a21 ) / det
  rsq = xc * xc + yc * yc
  diff = ( ( x0 - x1 - xc)**2 + ( y0 - y1 - yc )**2 ) - rsq
  tolabs = tol * rsq

  if ( diff < - tolabs ) then
    cmcirc = 1
  else if ( tolabs < diff ) then
    cmcirc = -1
  else
    cmcirc = 0
  end if

  return
end
subroutine cvdec2 ( angspc, angtol, nvc, npolg, nvert, maxvc, maxhv, &
  maxpv, maxiw, maxwk, vcl, regnum, hvl, pvl, iang, iwk, wk, ierror )

!*****************************************************************************80
!
!! CVDEC2 decomposes a polygonal region into convex polygons.
!
!  Purpose: 
!
!    Decompose general polygonal region (which is decomposed
!    into simple polygons on input) into convex polygons using
!    vertex coordinate list, head vertex list, and polygon vertex
!    list data structures.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ANGSPC, the angle spacing parameter in radians 
!    used in controlling vertices to be considered as an endpoint of a
!    separator.
!
!    Input, real ( kind = 8 ) ANGTOL, the angle tolerance parameter in radians 
!    used in accepting separator(s).
!
!    Input/output, integer NVC, the number of vertex coordinates or positions 
!    used in VCL.
!
!    Input/output, integer NPOLG, the number of polygonal subregions or 
!    positions used in HVL array.
!
!    Input/output, integer NVERT, the number of polygon vertices or positions 
!    used in PVL array.
!
!    Input, integer MAXVC, the maximum size available for VCL array, 
!    should be greater than or equal to the number of vertex coordinates
!    required for decomposition.
!
!    Input, integer MAXHV, the maximum size available for HVL, REGNUM arrays, 
!    should be greater than or equal to the number of polygons required for
!    decomposition.
!
!    Input, integer MAXPV, the maximum size available for PVL, IANG arrays; 
!    should be greater than or equal to the number of polygon vertices 
!    required for decomposition.
!
!    Input, integer MAXIW, the maximum size available for IWK array; should be 
!    about 3 times maximum number of vertices in any polygon.
!
!    Input, integer MAXWK, the maximum size available for WK array; should be 
!    about 5 times maximum number of vertices in any polygon.
!
!    Input/output, real ( kind = 8 ) VCL(1:2,1:NVC), the vertex coordinate list.
!
!    Input/output, integer REGNUM(1:NPOLG), region numbers.
!
!    Input/output, integer HVL(1:NPOLG), the head vertex list.
!
!    Input/output, integer PVL(1:4,1:NVERT), real ( kind = 8 ) IANG(1:NVERT),
!    the polygon vertex list and interior angles; see routine DSPGDC for 
!    more details.  Note that the data structures should be as output from 
!    routine SPDEC2.
!
!    Workspace, integer IWK(1:MAXIW).
!
!    Workspace, real ( kind = 8 ) WK(1:MAXWK).
!
!    Output, integer IERROR, error flag.  For abnormal return,
!    IERROR is set to 3, 4, 5, 6, 7, 206, 207, 208, 209, 210, or 212.
!
  integer maxhv
  integer maxiw
  integer maxpv
  integer maxvc
  integer maxwk

  real ( kind = 8 ) angspc
  real ( kind = 8 ) angtol
  integer hvl(maxhv)
  real ( kind = 8 ) iang(maxpv)
  integer ierror
  integer iwk(maxiw)
  integer npolg
  integer nvc
  integer nvert
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) piptol
  integer pvl(4,maxpv)
  integer regnum(maxhv)
  real ( kind = 8 ) tol
  integer v
  real ( kind = 8 ) vcl(2,maxvc)
  integer w1
  integer w2
  real ( kind = 8 ) wk(maxwk)

  ierror = 0
  tol = 100.0D+00 * epsilon ( tol )
!
!  For each reflex vertex, resolve it with one or two separators
!  and update VCL, HVL, PVL, IANG.
!
  piptol = pi + tol
  v = 1

  do

    if ( nvert < v ) then
      exit
    end if

    if ( piptol < iang(v) ) then

      call resvrt ( v, angspc, angtol, nvc, nvert, maxvc, maxpv, maxiw, &
        maxwk, vcl, pvl, iang, w1, w2, iwk, wk, ierror )

      if ( ierror /= 0 ) then
        return
      end if

      call insed2 ( v ,w1, npolg, nvert, maxhv, maxpv, vcl, regnum, hvl, &
        pvl, iang, ierror )

      if ( ierror /= 0 ) then
        return
      end if

      if ( 0 < w2 ) then
        call insed2 ( v, w2, npolg, nvert, maxhv, maxpv, vcl, regnum, hvl, &
          pvl, iang, ierror )
      end if

      if ( ierror /= 0 ) then
        return
      end if

    end if

    v = v + 1

  end do

  return
end
subroutine cvdtri ( inter, ldv, nt, vcl, til, tedg, sptr, ierror )

!*****************************************************************************80
!
!! CVDTRI converts boundary triangles to Delaunay triangles.
!
!  Purpose: 
!
!    Convert triangles in strip near boundary of polygon
!    or inside polygon to Delaunay triangles.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, logical INTER, is .TRUE. if and only if there is at least 
!    one interior mesh vertex.
!
!    Input, integer LDV, the leading dimension of VCL in calling routine.
!
!    Input, integer NT, the number of triangles in strip or polygon.
!
!    Input, VCL(1:2,1:*), the vertex coordinate list.
!
!    Input/output, integer TIL(1:3,1:NT), the triangle incidence list.
!
!    Input/output, integer TEDG(1:3,1:NT) - TEDG(J,I) refers to edge with 
!    vertices TIL(J:J+1,I) and contains index of merge edge or
!    a value greater than NT for edge of chains.
!
!    Workspace, SPTR(1:NT) - SPTR(I) = -1 if merge edge I is not in LOP stack,
!    else greater than or equal to 0 and pointer (index of SPTR) to next 
!    edge in stack (0 indicates bottom of stack).
!
!    Output, integer IERROR, error flag.  On abnormal return:
!    IERROR is set to 231.
!
  integer ldv
  integer nt

  integer e
  integer ierror
  integer ind(2)
  logical inter
  integer itr(2)
  integer k
  integer mxtr
  logical sflag
  integer sptr(nt)
  integer tedg(3,nt)
  integer til(3,nt)
  integer top
  real ( kind = 8 ) vcl(ldv,*)

  ierror = 0
  sflag = .true.
  sptr(1:nt) = -1

  do k = 1, nt

    mxtr = k + 1

    if ( k == nt ) then
      if ( .not. inter ) then
        return
      end if
      mxtr = nt
      sflag = .false.
    end if

    top = k
    sptr(k) = 0

    do

      e = top
      top = sptr(e)

      call fndtri ( e, mxtr, sflag, tedg, itr, ind, ierror )

      if ( ierror /= 0 ) then
        return
      end if

      call lop ( itr, ind, k, top, ldv, vcl, til, tedg, sptr )

      if ( top <= 0 ) then
        exit
      end if

    end do

  end do

  return
end
function degrees_to_radians ( angle )

!*****************************************************************************80
!
!! DEGREES_TO_RADIANS converts an angle from degrees to radians.
!
!  Modified:
!
!    10 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ANGLE, an angle in degrees.
!
!    Output, real ( kind = 8 ) DEGREES_TO_RADIANS, the equivalent angle
!    in radians.
!
  real ( kind = 8 ) angle
  real ( kind = 8 ) degrees_to_radians
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00

  degrees_to_radians = ( angle / 180.0D+00 ) * pi

  return
end
subroutine delaunay_print ( num_pts, xc, num_tri, nodtri, tnbr )

!*****************************************************************************80
!
!! DELAUNAY_PRINT prints out information defining a Delaunay triangulation.
!
!  Modified:
!
!    08 July 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NUM_PTS, the number of points.
!
!    Input, real ( kind = 8 ) XC(2,NUM_PTS), the point coordinates.
!
!    Input, integer NUM_TRI, the number of triangles.
!
!    Input, integer NODTRI(3,NUM_TRI), the nodes that make up the triangles.
!
!    Input, integer TNBR(3,NUM_TRI), the triangle neighbors on each side.
!
  integer num_pts
  integer num_tri

  integer i
  integer i4_wrap
  integer j
  integer k
  integer n1
  integer n2
  integer nodtri(3,num_tri)
  integer s
  integer t
  integer tnbr(3,num_tri)
  real ( kind = 8 ) xc(2,num_pts)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DELAUNAY_PRINT'
  write ( *, '(a)' ) '  Information defining a Delaunay triangulation.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The number of points is ', num_pts

  call r8mat_print ( num_pts, num_pts, 2, transpose ( xc ), &
    '  Point coordinates (transpose of internal array)' )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The number of triangles is ', num_tri
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Sets of three points are used as vertices of'
  write ( *, '(a)' ) '  the triangles.  For each triangle, the points'
  write ( *, '(a)' ) '  are listed in counterclockwise order.'

  call i4mat_print ( num_tri, num_tri, 3, transpose ( nodtri ), &
    '  Nodes that make up triangles (transpose of internal array)' )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  On each side of a given triangle, there is either'
  write ( *, '(a)' ) '  another triangle, or a piece of the convex hull.'
  write ( *, '(a)' ) '  For each triangle, we list the indices of the three'
  write ( *, '(a)' ) '  neighbors, or (if negative) the codes of the'
  write ( *, '(a)' ) '  segments of the convex hull.'

  call i4mat_print ( num_tri, num_tri, 3, transpose ( tnbr ), &
    '  Indices of neighboring triangles (transpose of internal array)' )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  The number of boundary points (and segments) is ', &
    2 * num_pts - num_tri - 2

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The segments that make up the convex hull can be'
  write ( *, '(a)' ) '  determined from the negative entries of the triangle'
  write ( *, '(a)' ) '  neighbor list.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  # Tri Side  N1  N2'
  write ( *, '(a)' ) ' '
  k = 0

  do i = 1, num_tri
    do j = 1, 3
      if ( tnbr(j,i) < 0 ) then
        s = - tnbr(j,i)
        t = s / 3
        s = mod ( s, 3 ) + 1
        k = k + 1
        n1 = nodtri(s,t)
        n2 = nodtri(i4_wrap(s+1,1,3),t)
        write ( *, '(5i4)' ) k, t, s, n1, n2
      end if
    end do
  end do

  return
end
subroutine dhpsrt ( k, n, lda, a, map )

!*****************************************************************************80
!
!! DHPSRT sorts points into lexicographic order using heap sort
!
!  Discussion:
!
!    This routine uses heapsort to obtain the permutation of N K-dimensional
!    points so that the points are in lexicographic increasing order.
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Modified:
!
!    19 February 2001
!
!  Parameters:
!
!    Input, integer K, the dimension of the points (for instance, 2
!    for points in the plane).
!
!    Input, integer N, the number of points.
!
!    Input, integer LDA, the leading dimension of array A in the calling
!    routine; LDA should be at least K.
!
!    Input, real ( kind = 8 ) A(LDA,N); A(I,J) contains the I-th coordinate
!    of point J.
!
!    Input/output, integer MAP(N).
!    On input, the points of A with indices MAP(1), MAP(2), ...,
!    MAP(N) are to be sorted.
!
!    On output, MAP contains a permutation of its input values, with the
!    property that, lexicographically,
!      A(*,MAP(1)) <= A(*,MAP(2)) <= ... <= A(*,MAP(N)).
!
  integer lda
  integer n

  real ( kind = 8 ) a(lda,*)
  integer i
  integer k
  integer map(n)

  do i = n/2, 1, -1
    call dsftdw ( i, n, k, lda, a, map )
  end do

  do i = n, 2, -1
    call i4_swap ( map(1), map(i) )
    call dsftdw ( 1, i-1, k, lda, a, map )
  end do

  return
end
function diaedg ( x0, y0, x1, y1, x2, y2, x3, y3 )

!*****************************************************************************80
!
!! DIAEDG triangulates 4 points using the circumcircle criterion.
!
!  Diagram:
!
!    3---2      3---2
!    |\  |      |  /|
!    | \ |  or  | / |
!    |  \|      |/  |
!    0---1      0---1
!
!  Discussion:
!
!    Given four points, to be organized into two triangles, the routine
!    chooses 0--2 or 1--3 as the diagonal edge, based on the circumcircle
!    criterion.  
!
!    The points may be regarded as the vertices of a simple quadrilateral,
!    and should be listed in counterclockwise order.
!
!  Modified:
!
!    13 April 2004
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X0, Y0, X1, Y1, X2, Y2, X3, Y3, the points.
!
!    Output, integer DIAEDG, chooses a diagonal:
!    +1, choose edge 0--2;
!    -1, choose edge 1--3;
!     0, the four vertices are essentially cocircular.
!
  implicit none

  real ( kind = 8 ) ca
  real ( kind = 8 ) cb
  integer diaedg
  real ( kind = 8 ) dx10
  real ( kind = 8 ) dx12
  real ( kind = 8 ) dx30
  real ( kind = 8 ) dx32
  real ( kind = 8 ) dy10
  real ( kind = 8 ) dy12
  real ( kind = 8 ) dy30
  real ( kind = 8 ) dy32
  real ( kind = 8 ) s
  real ( kind = 8 ) tol
  real ( kind = 8 ) tola
  real ( kind = 8 ) tolb
  real ( kind = 8 ) x0
  real ( kind = 8 ) x1
  real ( kind = 8 ) x2
  real ( kind = 8 ) x3
  real ( kind = 8 ) y0
  real ( kind = 8 ) y1
  real ( kind = 8 ) y2
  real ( kind = 8 ) y3

  tol = 100.0D+00 * epsilon ( tol )

  dx10 = x1 - x0
  dy10 = y1 - y0
  dx12 = x1 - x2
  dy12 = y1 - y2
  dx30 = x3 - x0
  dy30 = y3 - y0
  dx32 = x3 - x2
  dy32 = y3 - y2

  tola = tol * max ( abs ( dx10 ), abs ( dy10 ), abs ( dx30 ), abs ( dy30 ) )
  tolb = tol * max ( abs ( dx12 ), abs ( dy12 ), abs ( dx32 ), abs ( dy32 ) )

  ca = dx10 * dx30 + dy10 * dy30
  cb = dx12 * dx32 + dy12 * dy32
!
!  If both angles 301 and 123 are acute, choose 1-3 as the edge.
!
  if ( tola < ca .and. tolb < cb ) then

    diaedg = -1
!
!  If both angles 301 and 123 are obtuse, choose 0-2 as the edge.
!
  else if ( ca < -tola .and. cb < -tolb ) then

    diaedg = 1
 
  else

    tola = max ( tola, tolb )
    s = ( dx10 * dy30 - dx30 * dy10 ) * cb + ( dx32 * dy12 - dx12 * dy32 ) * ca

    if ( tola < s ) then
      diaedg = -1
    else if ( s < -tola ) then
      diaedg = 1
    else
      diaedg = 0
    end if

  end if

  return
end
subroutine diam2 ( nvrt, xc, yc, i1, i2, diamsq, ierror )

!*****************************************************************************80
!
!! DIAM2 finds the diameter of a convex polygon.
!
!  Purpose: 
!
!    Find the diameter of a convex polygon with vertices
!    given in counter clockwise order and with all interior angles < PI.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer NVRT, the number of vertices.
!
!    Input, real ( kind = 8 ) XC(NVRT),YC(NVRT), the vertex coordinates in 
!    counter-clockwise order.
!
!    Output, integer I1, I2 , indices in XC, YC of the diameter edge; the
!    diameter is from (XC(I1),YC(I1)) to (XC(I2),YC(I2)).
!
!    Output, real ( kind = 8 ) DIAMSQ, the square of the diameter.
!
!    Output, integer IERROR, an error flag.
!    0, no error was detected.
!    200, an error was detected.
!
  integer nvrt

  real ( kind = 8 ) area1
  real ( kind = 8 ) area2
  real ( kind = 8 ) areatr
  real ( kind = 8 ) c1mtol
  real ( kind = 8 ) c1ptol
  real ( kind = 8 ) diamsq
  real ( kind = 8 ) dist
  integer i1
  integer i2
  integer ierror
  integer j
  integer jp1
  integer k
  integer kp1
  integer m
  real ( kind = 8 ) tol
  real ( kind = 8 ) xc(nvrt)
  real ( kind = 8 ) yc(nvrt)
!
  ierror = 0
  tol = 100.0D+00 * epsilon ( tol )
!
!  Find the first vertex which is farthest from the edge connecting
!  vertices with indices NVRT, 1.
!
  c1mtol = 1.0D+00 - tol
  c1ptol = 1.0D+00 + tol
  j = nvrt
  jp1 = 1
  k = 2
  area1 = areatr ( xc(j), yc(j), xc(jp1), yc(jp1), xc(k), yc(k) )

  do

    area2 = areatr ( xc(j), yc(j), xc(jp1), yc(jp1), xc(k+1), yc(k+1) )

    if ( area2 <= area1 * c1ptol ) then
      exit
    end if

    area1 = area2
    k = k + 1

  end do

  m = k
  diamsq = 0.0D+00
!
!  Find diameter = maximum distance of antipodal pairs.
!
  area1 = areatr ( xc(j), yc(j), xc(jp1), yc(jp1), xc(k), yc(k) )

  do

    kp1 = k + 1
    if ( nvrt < kp1 ) then
      kp1 = 1
    end if

    area2 = areatr ( xc(j), yc(j), xc(jp1), yc(jp1), xc(kp1), yc(kp1) )

    if ( area1 * c1ptol < area2 ) then
      k = k + 1
      area1 = area2
    else if ( area2 < area1 * c1mtol ) then
      j = jp1
      jp1 = j + 1
      area1 = areatr ( xc(j), yc(j), xc(jp1), yc(jp1), xc(k), yc(k) )
    else
      k = k + 1
      j = jp1
      jp1 = j + 1
      area1 = areatr ( xc(j), yc(j), xc(jp1), yc(jp1), xc(k), yc(k) )
    end if

    if ( m < j .or. nvrt < k ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DIAM2 - Fatal error!'
      write ( *, '(a)' ) '  M < J or NVRT < K.'
      ierror = 200
      return
    end if

    dist = ( xc(j) - xc(k) )**2 + ( yc(j) - yc(k) )**2

    if ( diamsq < dist ) then
      diamsq = dist
      i1 = j
      i2 = k
    end if

    if ( j == m .and. k == nvrt ) then
      exit
    end if

  end do

  return
end
function dless ( k, p, q )

!*****************************************************************************80
!
!! DLESS determine whether P is lexicographically less than Q.
!
!  Discussion:
!
!    P and Q are K-dimensional points.
!
!  Modified:
!
!    19 February 2001
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer K, the spatial dimension of the points.
!
!    Input, real ( kind = 8 ) P(K), Q(K), the points to be compared.
!
!    Output, logical RLESS, is TRUE if P < Q, FALSE otherwise.
!
  integer k

  real ( kind = 8 ) cmax
  logical dless
  integer i
  real ( kind = 8 ) p(k)
  real ( kind = 8 ) q(k)
  real ( kind = 8 ) tol

  tol = 100.0D+00 * epsilon ( tol )

  do i = 1, k

    cmax = max ( abs ( p(i) ), abs ( q(i) ) )

    if ( tol * cmax < abs ( p(i) - q(i) ) .and. tol < cmax ) then

      if ( p(i) < q(i) ) then
        dless = .true.
      else
        dless = .false.
      end if

      return
    end if

  end do

  dless = .false.

  return
end
subroutine dsftdw ( l, u, k, lda, a, map )

!*****************************************************************************80
!
!! DSFTDW sifts A(*,MAP(L)) down a heap of size U.
!
!  Modified:
!
!    19 February 2001
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer L, U, the lower and upper indices of part of the heap.
!
!    Input, integer K, the spatial dimension of the points.
!
!    Input, integer LDA, the leading dimension of A in the calling routine.
!
!    Input, real ( kind = 8 ) A(LDA,N); A(I,J) contains the I-th coordinate
!    of point J.
!
!    Input/output, integer MAP(N).
!    On input, the points of A with indices MAP(1), MAP(2), ...,
!    MAP(N) are to be sorted.
!
!    On output, MAP contains a permutation of its input values, with the
!    property that, lexicographically,
!      A(*,MAP(1)) <= A(*,MAP(2)) <= ... <= A(*,MAP(N)).
!
  integer lda

  real ( kind = 8 ) a(lda,*)
  logical dless
  integer i
  integer j
  integer k
  integer l
  integer map(*)
  integer t
  integer u

  i = l
  j = 2 * i
  t = map(i)

  do while ( j <= u )

    if ( j < u ) then
      if ( dless ( k, a(1,map(j)), a(1,map(j+1)) ) ) then
        j = j + 1
      end if
    end if

    if ( dless ( k, a(1,map(j)), a(1,t) ) ) then
      exit
    end if

    map(i) = map(j)
    i = j
    j = 2 * i

  end do

  map(i) = t

  return
end
subroutine dsmcpr ( nhole, nvbc, vcl, maxhv, maxpv, maxho, nvc, npolg, &
  nvert, nhola, regnum, hvl, pvl, iang, holv, ierror )

!*****************************************************************************80
!
!! DSMCPR initializes the polygonal decomposition data structure.
!
!  Purpose: 
!
!    Initialize the polygonal decomposition data structure
!    given a multiply-connected polygonal region with 1 outer
!    boundary curve and 0 or more inner boundary curves of holes.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer NHOLE, the number of holes in the region.
!
!    Input, integer NVBC(1:NHOLE+1), the number of vertices per boundary curve; 
!    first boundary curve is the outer boundary of the region.
!
!    Input, real ( kind = 8 ) VCL(1:2,1:NVC), the vertex coordinates of 
!    boundary curves in counter clockwise order; 
!    NVC = NVBC(1) + ... + NVBC(NHOLE+1); 
!    positions 1 to NVBC(1) of VCL contain the vertex coordinates of the
!    outer boundary in counter clockwise order; positions NVBC(1)+1 to
!    NVBC(1)+NVBC(2) contain the vertex coordinates of the
!    first hole boundary in counter clockwise order, etc.
!
!    Input, integer MAXHV, the maximum size available for HVL, REGNUM arrays, 
!    should be greater than or equal to NHOLE + 1.
!
!    Input, integer MAXPV, the maximum size available for PVL, IANG arrays;
!    should be greater than or equal to NVC.
!
!    Input, integer MAXHO, the maximum size available for HOLV array; should be
!    greater than or equal to NHOLE*2.
!
!    Output, integer NVC, the number of vertex coordinates, set to sum of 
!    NVBC(I).
!
!    Output, integer NPOLG, the number of polygonal subregions, set to 1.
!    For consistency with DSPGDC.
!
!    Output, integer NVERT, the number of vertices in PVL, set to NVC.
!    For consistency with DSPGDC.
!
!    Output, integer NHOLA, number of attached holes, set to 0.
!    For consistency with DSPGDC.
!
!    Output, integer REGNUM(1:1), region number of only subregion, set to 1
!    For consistency with DSPGDC.
!
!    Output, integer HVL(1:NHOLE+1), the head vertex list; first entry is the 
!    head vertex (index in PVL) of outer boundary curve; next
!    NHOLE entries contain the head vertex of a hole.
!
!    Output, integer PVL(1:4,1:NVC), IANG(1:NVC), the polygon vertex list and 
!    interior angles; vertices of outer boundary curve are in counter clockwise 
!    order followed by vertices of each hole in CW hole; vertices
!    of each polygon are in a circular linked list; see
!    routine DSPGDC for more details of this data structure.
!
!    Output, integer HOLV(1:NHOLE*2), the indices in PVL of top and bottom 
!    vertices of holes; first (last) NHOLE entries are for top (bottom)
!    vertices; top (bottom) vertices are sorted in decreasing
!    (increasing) lexicographic (y,x) order of coordinates.
!
!    Output, integer IERROR, error flag.
!    For abnormal return, IERROR is set to 2, 4, or 5.
!
  integer maxho
  integer maxpv
  integer nhole

  real ( kind = 8 ) angle
  integer, parameter :: edgv = 4
  integer i
  real ( kind = 8 ) iang(maxpv)
  integer ierror
  integer iv
  integer ivs
  integer j
  integer, parameter :: loc = 1
  integer lv
  integer lvp
  integer lvs
  integer maxhv
  integer nhola
  integer npolg
  integer nv
  integer nvc
  integer nvert
  integer nvs
  integer hvl(nhole+1)
  integer holv(maxho)
  integer nvbc(nhole+1)
  integer, parameter :: polg = 2
  integer pvl(4,maxpv)
  integer regnum(1)
  integer, parameter :: succ = 3
  real ( kind = 8 ) vcl(2,*)

  ierror = 0
  nvc = sum ( nvbc(1:nhole+1) )
  npolg = 1
  nvert = nvc
  nhola = 0
  regnum(1) = 1

  if ( maxhv < nhole + 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DSMCPR - Fatal error!'
    write ( *, '(a)' ) '  MAXHV < NHOLE + 1.'
    ierror = 4
    return
  end if

  if ( maxpv < nvc ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DSMCPR - Fatal error!'
    write ( *, '(a)' ) '  MAXPV < NVC.'
    ierror = 5
    return
  end if

  if ( maxho < 2 * nhole ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DSMCPR - Fatal error!'
    write ( *, '(a)' ) '  MAXHO < 2 * NHOLE.'
    ierror = 2
    return
  end if
!
!  Initialize the HVL and PVL arrays.
!
  hvl(1) = 1
  nv = nvbc(1)

  do i = 1, nv
    pvl(loc,i) = i
    pvl(polg,i) = 1
    pvl(succ,i) = i + 1
    pvl(edgv,i) = 0
  end do

  pvl(succ,nv) = 1

  do j = 1, nhole
    hvl(j+1) = nv + 1
    nvs = nv + nvbc(j+1)
    do i = nv+1, nvs
      pvl(loc,i) = i
      pvl(polg,i) = 1
      pvl(succ,i) = i - 1
      pvl(edgv,i) = 0
    end do
    pvl(succ,nv+1) = nvs
    nv = nvs
  end do
!
!  Initialize the IANG array.
!
  do i = 1, nhole+1

    j = hvl(i)
    lvp = pvl(loc,j)
    iv = pvl(succ,j)
    lv = pvl(loc,iv)

    do

      ivs = pvl(succ,iv)
      lvs = pvl(loc,ivs)
      iang(iv) = angle ( vcl(1,lvp), vcl(2,lvp), vcl(1,lv), vcl(2,lv), &
        vcl(1,lvs), vcl(2,lvs) )

      if ( iv == j ) then
        exit
      end if

      lvp = lv
      iv = ivs
      lv = lvs

    end do

  end do
!
!  Initialize the HOLV array.
!
  if ( 0 < nhole ) then
    call holvrt ( nhole, vcl, hvl(2), pvl, holv )
  end if

  return
end
subroutine dsmdf2 ( hflag, nvc, npolg, maxwk, vcl, hvl, pvl, iang, ivrt, &
  xivrt, widsq, edgval, vrtval, area, wk, ierror )

!*****************************************************************************80
!
!! DSMDF2 sets up a data structure for a heuristic mesh distribution.
!
!  Purpose: 
!
!    Set up the data structure for heuristic mesh distribution
!    function from data structure for convex polygon decomposition
!    if HFLAG is .TRUE., else set up only IVRT and XIVRT.
!
!    Also compute areas of convex polygons.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, logical HFLAG, set to .TRUE. if data structure is to be constructed,
!    .FALSE. if only IVRT, XIVRT, AREA are to be computed.
!
!    Input, integer NVC, the number of vertex coordinates in VCL array.
!
!    Input, integer NPOLG, the number of polygonal subregions in HVL array.
!
!    Input, integer MAXWK, the maximum size available for WK array; should be
!    2 times maximum number of vertices in any polygon.
!
!    Input, VCL(1:2,1:NVC), the vertex coordinate list.
!
!    Input, integer HVL(1:NPOLG), the head vertex list.
!
!    Input, integer PVL(1:4,1:*), real ( kind = 8 ) IANG(1:*), the polygon
!    vertex list, interior angles.
!
!    Output, integer IVRT(1:*), the indices of polygon vertices in VCL, ordered 
!    by polygon; same size as PVL.  For heuristic MDF data structure.
!
!    Output, XIVRT(1:NPOLG+1), the pointer to first vertex of each polygon
!    in IVRT; vertices of polygon K are IVRT(I) for I from
!    XIVRT(K) to XIVRT(K+1)-1.  For heuristic MDF data structure.
!
!    Output, real ( kind = 8 ) WIDSQ(1:NPOLG), the square of width of convex 
!    polygons.  For heuristic MDF data structure.
!
!    Output, real ( kind = 8 ) EDGVAL(1:*), the value associated with each 
!    edge of decomposition; same size as PVL.  For heuristic MDF data structure.
!
!    Output, real ( kind = 8 ) VRTVAL(1:NVC), the value associated with each 
!    vertex of decomposition.  For heuristic MDF data structure.
!
!    Output, real ( kind = 8 ) AREA(1:NPOLG), the area of convex polygons.
!
!    Workspace, real ( kind = 8 ) WK(1:MAXWK).
!
!    Output, integer IERROR, error flag.  On abnormal return,
!    IERROR is set to 7 or 201.
!
  integer maxwk
  integer npolg
  integer nvc

  real ( kind = 8 ) area(npolg)
  real ( kind = 8 ) areapg
  integer, parameter :: edgv = 4
  real ( kind = 8 ) edgval(*)
  logical hflag
  integer hvl(npolg)
  integer i
  real ( kind = 8 ) iang(*)
  integer ierror
  integer il
  integer ivrt(*)
  integer j
  integer jl
  integer k
  integer l
  integer, parameter :: loc = 1
  integer m
  integer nvrt
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) pimtol
  integer, parameter :: polg = 2
  integer pvl(4,*)
  real ( kind = 8 ) s
  integer, parameter :: succ = 3
  real ( kind = 8 ) tol
  real ( kind = 8 ) vcl(2,nvc)
  real ( kind = 8 ) vrtval(nvc)
  real ( kind = 8 ) widsq(npolg)
  real ( kind = 8 ) wk(maxwk)
  integer xc
  integer xivrt(npolg+1)
  integer yc

  ierror = 0
  tol = 100.0D+00 * epsilon ( tol )
!
!  Compute area and square of width of polygons.
!
  pimtol = pi - tol

  do k = 1, npolg

    nvrt = 0
    i = hvl(k)

    do

      if ( iang(i) < pimtol ) then
        nvrt = nvrt + 1
      end if

      i = pvl(succ,i)

      if ( i == hvl(k) ) then
        exit
      end if

    end do

    if ( maxwk < 2 * nvrt ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DSMDF2 - Fatal error!'
      write ( *, '(a)' ) '  MAXWK < 2 * NVRT.'
      write ( *, '(a,i6)' ) '  NVRT = ', nvrt
      write ( *, '(a,i6)' ) '  MAXWK = ', maxwk
      ierror = 7
      return
    end if

    xc = 0

    do

      if ( iang(i) < pimtol ) then
        j = pvl(loc,i)
        xc = xc + 1
        wk(xc) = vcl(1,j)
        wk(xc+nvrt) = vcl(2,j)
      end if

      i = pvl(succ,i)

      if ( i == hvl(k) ) then
        exit
      end if

    end do

    xc = 1
    yc = xc + nvrt
    area(k) = areapg ( nvrt, wk(xc), wk(yc) ) * 0.5D+00

    if ( hflag ) then

      call width2 ( nvrt, wk(xc), wk(yc), i, j, widsq(k), ierror )

      if ( ierror /= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'DSMDF2 - Fatal error!'
        write ( *, '(a)' ) '  WIDTH2 returns an error condition.'
        return
      end if

    end if

  end do
!
!  Set up IVRT, XIVRT, EDGVAL, VRTVAL arrays.
!
  l = 1

  do k = 1, npolg

    xivrt(k) = l
    i = hvl(k)
    il = pvl(loc,i)

    do

      ivrt(l) = il
      j = pvl(succ,i)
      jl = pvl(loc,j)

      if ( hflag ) then

        s = min ( (vcl(1,jl) - vcl(1,il) )**2 + ( vcl(2,jl) - vcl(2,il) )**2, &
          widsq(k) )

        m = pvl(edgv,i)
        if ( 0 < m ) then
          s = min ( s, widsq(pvl(polg,m) ) )
        end if

        edgval(l) = s

      end if

      l = l + 1
      i = j
      il = jl

      if ( i == hvl(k) ) then
        exit
      end if

    end do

  end do

  xivrt(npolg+1) = l

  if ( .not. hflag ) then
    return
  end if

  vrtval(1:nvc) = 0.0D+00

  do k = 1, npolg

    j = xivrt(k+1) - 1
    l = j

    do i = xivrt(k),l

      il = ivrt(i)

      if ( vrtval(il) == 0.0D+00 ) then
        vrtval(il) = min ( edgval(i), edgval(j) )
      else
        vrtval(il) = min ( vrtval(il), edgval(i), edgval(j) )
      end if

      j = i

    end do

  end do

  return
end
subroutine dspgdc ( nvc, vcl, incr, ncur, nvbc, icur, ivrt, maxhv, maxpv,  &
  maxho, npolg, nvert, nhole, nhola, regnum, hvl, pvl, iang, holv, htsiz,  &
  maxedg, ht, edge, map, ierror )

!*****************************************************************************80
!
!! DSPGDC initializes the polygonal decomposition data structure.
!
!  Purpose: 
!
!    Initialize the polygonal decomposition data structure
!    given an initial decomposition of a polygonal region which
!    may have holes and/or cut, separator, and hole interfaces.
!    Holes and hole interfaces must be simple polygons.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer NVC, the number of distinct vertex coordinates in region.
!
!    Input, real ( kind = 8 ) VCL(1:2,1:NVC), the vertex coordinates of 
!    boundary curves in arbitrary order.
!
!    Input, integer INCR, a positive integer greater than or equal to NVC, 
!    e.g. 10000, added to some elements of IVRT array.
!
!    Input, integer NCUR, the number of boundary curves (includes outer boundary
!    curves of subregions and boundary curves of holes
!    and hole interfaces).
!
!    Input, integer NVBC(1:NCUR), the number of vertices per boundary curve.
!
!    Input, integer ICUR(1:NCUR), indicates type and location of the curves:
!    ICUR(I) = 0 if Ith curve is outer boundary curve,
!    ICUR(I) = K if Ith curve is a hole and is inside
!      the subregion to the left of Kth curve,
!    ICUR(I) = -K if Ith curve is a hole interface and is
!      inside the subregion to the left of Kth curve.
!    K must be the index of an outer or hole interface
!    boundary curve (hole interfaces may be nested).
!    If the Ith curve is inside more than one subregion
!    due to nesting of hole interfaces, then the subregion
!    to the left of Kth curve must be the smallest
!    subregion containing the Ith curve.
!
!    Input, integer IVRT(1:NV), indices in VCL of vertices of boundary curves;
!    NV = NVBC(1) + ... + NVBC(NCUR); the vertices of each
!    boundary curve must be in counter clockwise order; the first NVBC(1)
!    positions of IVRT are used for the first curve; the
!    next NVBC(2) positions are used for second curve, etc.
!    If the Ith curve is the outer boundary of a subregion
!    determined from cut and separator interfaces, then the
!    elements of IVRT which correspond to this curve are used
!    both for an index in VCL and indicating the type of the
!    edge joining a vertex and its successor as follows.
!    Let J be in range of positions used for the Ith curve
!    and K be the index in VCL of the coordinates of a vertex
!    of the Ith curve. Consider the edge originating from this
!    vertex. IVRT(J) = -K if the edge is part of a cut or
!    separator interface (i.e. there is a subregion to right
!    of edge). IVRT(J) = K if the edge is part of the outer
!    boundary of the region (i.e. the unbounded exterior of
!    the region is to the right of edge). IVRT(J) = K + INCR
!    if the edge is part of the boundary of a hole (i.e.
!    there is a bounded area to the right of edge which is
!    not in the region. If the Ith curve is the boundary of
!    a hole or hole interface, then only IVRT(J) = K is used.
!
!    Input, integer MAXHV, the maximum size available for HVL, REGNUM arrays, 
!    should be greater than or equal to NCUR + (number of hole interfaces).
!
!    Input, integer MAXPV, the maximum size available for PVL, IANG arrays; 
!    should be greater than or equal to NVERT (see below).
!
!    Input, integer MAXHO, the maximum size available for HOLV array; should be
!    greater than or equal to NHOLE*2 + NHOLA (see below).
!
!    Input, integer HTSIZ, the size of hash table HT; should be a prime number 
!    which is about NSC/2 where NSC is number of separator and cut
!    interface edges.
!
!    Input, integer MAXEDG, the maximum size available for EDGE array; should 
!    be at least NSC.
!
!    Output, integer NPOLG, the number of polygonal subregions, set to number 
!    of outer subregion boundaries plus number of hole interfaces.
!
!    Output, integer NVERT, the number of vertices in PVL, set to NV plus number
!    of vertices in holes and hole interfaces (< 2*NV).
!
!    Output, integer NHOLE, the number of holes and hole interfaces.
!
!    Output, integer NHOLA, the number of 'attached' holes; these holes are 
!    attached to the outer boundary of a subregion through vertices
!    or cut interfaces and have their edges in consecutive
!    order on the boundary (<= NV/4).
!
!    Output, integer REGNUM(1:NPOLG). region numbers to left of outer and hole
!    interface boundary curves, which are set to the indices
!    of ICUR or NVBC; this array may be useful in some
!    applications for identifying which original region a
!    subpolygon belongs to.
!
!    Output, HVL(1:NPOLG+NHOLE), the head vertex list; the first NPOLG
!    positions contain the head vertex (index in PVL) of an
!    outer or hole interface boundary curve in which the
!    vertices of the curve are in counter clockwise order in PVL; next
!    NHOLE positions contain the head vertex of a hole or
!    hole interface in which vertices are in CW order in PVL.
!
!    Output, PVL(1:4,1:NVERT), real ( kind = 8 ) IANG(1:NVERT), the polygon 
!    vertex list and interior angles; contains the 5 'arrays' LOC, POLG, SUCC
!    EDGV, IANG (the first 4 are integer arrays, the last
!    is a double precision array); the vertices of each
!    polygon (except for holes) are stored in counter clockwise order in a
!    circular linked list. PVL(LOC,V) is the location in VCL
!    of the coordinates of 'vertex' (index) V. IANG(V) is
!    the interior angle at vertex V. PVL(POLG,V) is polygon
!    number (index of HVL) of subregion containing vertex V
!    (this entry is different from the polygon index only
!    for holes). PVL(SUCC,V) is index in PVL of successor
!    vertex of vertex V. PVL(EDGV,V) gives information about
!    the edge joining vertices V and its successor - if the
!    edge is part of 1 polygon then PVL(EDGV,V) = 0; if the
!    edge is common to 2 polygons then 0 < PVL(EDGV,V) and
!    is equal to the index in PVL of the successor vertex
!    as represented in the other polygon; i.e. in latter
!    case, PVL(LOC,PVL(EDGV,V)) = PVL(LOC,PVL(SUCC,V)) and
!    PVL(EDGV,PVL(EDGV,V)) = V.
!
!    Output, integer HOLV(1:NHOLE*2+NHOLA), indices in PVL of top or bottom 
!    vertex of holes; first (next) NHOLE entries are for top (bottom)
!    vertices of holes and hole interfaces, with top (bottom)
!    vertices sorted in decreasing (increasing) lexicographic
!    (y,x) order of coord; last NHOLA entries are for attached
!    holes; if bottom vertex of attached hole is a simple
!    vertex of boundary curve containing the hole then entry
!    contains index of bottom vertex otherwise entry contains
!    index of top vertex (which is simple).
!
!    Workspace, integer MAP(1:NCUR), used for mapping input boundary curve 
!    numbers to polygon numbers.
!
!    Workspace, HT(0:HTSIZ-1), EDGE(1:4,1:MAXEDG) - hash table and edge records
!    used to determine matching occurrences of separator or
!    cut interface edges by calling routine EDGHT.
!
!    Output, integer IERROR, error flag.  On abnormal return,
!    IERROR is set to 1, 2, 4, 5, 215, or 216.
!
  integer htsiz
  integer maxedg
  integer maxho
  integer maxhv
  integer maxpv
  integer ncur

  real ( kind = 8 ) angle
  integer edge(4,maxedg)
  integer, parameter :: edgv = 4
  logical first
  integer hdfree
  integer holv(maxho)
  integer ht(0:htsiz-1)
  integer hvl(maxhv)
  integer i
  real ( kind = 8 ) iang(maxpv)
  integer icur(ncur)
  integer ierror
  integer incr
  integer ipoly
  integer iv
  integer ivrt(*)
  integer ivs
  integer j
  integer j1
  integer j2
  integer jend
  integer jstr
  integer k
  integer kmax
  integer kmin
  integer kpoly
  integer l
  integer last
  integer, parameter :: loc = 1
  integer lv
  integer lvp
  integer lvs
  integer map(ncur)
  integer mpoly
  integer nh2
  integer nhola
  integer nhole
  integer nholi
  integer nht
  integer npolg
  integer nv
  integer nvbc(ncur)
  integer nvc
  integer nvert
  integer, parameter :: polg = 2
  integer pvl(4,maxpv)
  integer regnum(maxhv)
  integer, parameter :: succ = 3
  real ( kind = 8 ) vcl(2,nvc)
  real ( kind = 8 ) x
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) y
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin

  ierror = 0
  nhola = 0
  nhole = 0
  nholi = 0
  nvert = 0

  do i = 1, ncur

    nvert = nvert + nvbc(i)

    if ( 0 < icur(i) ) then
      nhole = nhole + 1
    else if ( icur(i) < 0 ) then
      nholi = nholi + 1
      nvert = nvert + nvbc(i)
    end if

  end do

  npolg = ncur - nhole
  ipoly = 0
  iv = 0
  nv = 0
  hdfree = 0
  last = 0
  nht = 0

  ht(0:htsiz-1) = 0

  if ( maxhv < ncur + nholi ) then
    ierror = 4
    return
  else if ( maxpv < nvert ) then
    ierror = 5
    return
  else if ( maxho < ( nhole + nholi ) * 2 ) then
    ierror = 2
    return
  end if
!
!  Initialize REGNUM, HVL, PVL arrays for outer boundary curves.
!
  do i = 1, ncur

    if ( icur(i) /= 0 ) then
      map(i) = 0
      nv = nv + nvbc(i)
      cycle
    end if

    ipoly = ipoly + 1
    regnum(ipoly) = i
    hvl(ipoly) = iv + 1
    map(i) = ipoly
    jstr = nv + 1
    jend = nv + nvbc(i)

    do j = jstr, jend

      iv = iv + 1
      pvl(loc,iv) = abs ( ivrt(j) )
      pvl(polg,iv) = ipoly
      pvl(succ,iv) = iv + 1

      if ( 0 < ivrt(j) ) then
        pvl(edgv,iv) = 0
      else
!
!  The edge originating from current vertex is on a cut or
!  separator interface. Search in hash table for edge, and
!  insert or delete edge.  Set EDGV value if possible.
!
         lv = abs ( ivrt(j) )
         if ( incr < lv ) then
           lv = lv - incr
         end if

         if ( j < jend ) then
           lvs = abs ( ivrt(j+1) )
         else
           lvs = abs ( ivrt(jstr) )
         end if

         if ( incr < lvs ) then
           lvs = lvs - incr
         end if

         call edght ( lv, lvs, iv, nvc, htsiz, maxedg, hdfree, last, ht, &
           edge, ivs, ierror )

         if ( ierror /= 0 ) then
           return
         end if

         if ( 0 < ivs ) then
           pvl(edgv,iv) = ivs
           pvl(edgv,ivs) = iv
           nht = nht - 1
         else
           nht = nht + 1
         end if

       end if

     end do

     pvl(succ,iv) = hvl(ipoly)

     nv = nv + nvbc(i)

  end do

  if ( nht /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DSPGDC - Fatal error!'
    write ( *, '(a)' ) '  NHT /= 0.'
    ierror = 215
    return
  end if
!
!  Initialize REGNUM, HVL, PVL arrays for the hole interfaces.
!
  if ( nholi == 0 ) then
    go to 100
  end if

  do i = 1, ncur

    if ( icur(i) < 0 ) then
      ipoly = ipoly + 1
      map(i) = ipoly
    end if

  end do

  nv = 0

  do i = 1, ncur

    if ( icur(i) < 0 ) then

      ipoly = ipoly + 1
      kpoly = ipoly - nholi
      mpoly = map(-icur(i))
      regnum(kpoly) = i
      hvl(kpoly) = iv + 1
      hvl(ipoly) = iv + 2
      jstr = nv + 1
      jend = nv + nvbc(i)

      do j = jstr, jend
        iv = iv + 2
        pvl(loc,iv-1) = ivrt(j)
        pvl(polg,iv-1) = kpoly
        pvl(succ,iv-1) = iv + 1
        pvl(edgv,iv-1) = iv + 2
        pvl(loc,iv) = ivrt(j)
        pvl(polg,iv) = mpoly
        pvl(succ,iv) = iv - 2
        pvl(edgv,iv) = iv - 3
      end do

      pvl(succ,iv-1) = hvl(kpoly)
      pvl(edgv,iv-1) = hvl(ipoly)
      pvl(succ,hvl(ipoly)) = iv
      pvl(edgv,hvl(ipoly)) = iv - 1

    end if

    nv = nv + nvbc(i)
   
  end do
!
!  Initialize HVL, PVL arrays for the ordinary holes.
!
100 continue

  if ( nhole /= 0 ) then

    nv = 0

    do i = 1, ncur

      if ( 0 < icur(i) ) then

        ipoly = ipoly + 1
        mpoly = map(icur(i))
        hvl(ipoly) = iv + 1
        jstr = nv + 1
        jend = nv + nvbc(i)

        do j = jstr, jend
          iv = iv + 1
          pvl(loc,iv) = ivrt(j)
          pvl(polg,iv) = mpoly
          pvl(succ,iv) = iv - 1
          pvl(edgv,iv) = 0
        end do

        pvl(succ,hvl(ipoly)) = iv

      end if

      nv = nv + nvbc(i)

    end do

  end if
!
!  Determine bottom or top simple vertex of attached holes.
!
  nhole = nhole + nholi
  nh2 = nhole + nhole
  j1 = 0
  j2 = 0

  do i = 1, npolg-nholi

   j = hvl(i)

  150    continue

      if ( incr < pvl(loc,j) ) then
         j = pvl(succ,j)
         if ( j /= hvl(i) ) then
          go to 150
         else
          ierror = 216
          return
         end if
      end if

   first = .true.

  160    continue

      lv = pvl(loc,j)

      if ( 0 < j1 ) then
        if ( lv <= incr ) then
          j2 = j
        else if ( lv - incr == lvs ) then
          j2 = j
        else
          pvl(loc,j) = lv - incr
        end if
      else if ( incr < lv ) then
        j1 = j
        lvs = lv - incr
        pvl(loc,j) = lvs
      end if

      if ( 0 < j2 ) then
!
!  (Part of) hole starts at vertex J1 and ends at J2.
!
         if ( lv <= incr .and. lv /= lvs ) then
           go to 180
         end if

         k = j1

  170          continue

          if ( k == j1 ) then
            kmin = k
            kmax = k
            xmin = vcl(1,lvs)
            ymin = vcl(2,lvs)
            xmax = xmin
            ymax = ymin

          else

            l = pvl(loc,k)
            x = vcl(1,l)
            y = vcl(2,l)

            if ( y < ymin .or. y == ymin .and. x < xmin ) then
              kmin = k
              xmin = x
              ymin = y
            else if ( ymax < y .or. y == ymax .and. xmax < x ) then
              kmax = k
              xmax = x
              ymax = y
            end if

          end if

          k = pvl(succ,k)
         if ( k /= j2 ) then
           go to 170
         end if

         if ( kmin == j1 ) then
           kmin = kmax
         end if

         nhola = nhola + 1

         if ( maxho < nh2 + nhola ) then
          ierror = 2
          return
         end if

         holv(nh2+nhola) = kmin

180      continue

         j1 = 0
         j2 = 0

         if ( incr < lv ) then
           j1 = j
           pvl(loc,j) = lvs
         end if

      end if
      j = pvl(succ,j)

   if ( first ) then
     first = .false.
     jend = j
     go to 160
   else if ( j /= jend ) then
     go to 160
   end if

  end do
!
!  Initialize the IANG array.
!
  do i = 1, npolg+nhole

    j = hvl(i)
    lvp = pvl(loc,j)
    iv = pvl(succ,j)
    lv = pvl(loc,iv)

    do

      ivs = pvl(succ,iv)
      lvs = pvl(loc,ivs)
      iang(iv) = angle ( vcl(1,lvp), vcl(2,lvp), vcl(1,lv), vcl(2,lv), &
        vcl(1,lvs), vcl(2,lvs) )

      if ( iv == j ) then
        exit
      end if

      lvp = lv
      iv = ivs
      lv = lvs

    end do

  end do
!
!  Initialize HOLV array.
!
  if ( 0 < nhole ) then
    call holvrt ( nhole, vcl, hvl(npolg+1), pvl, holv )
  end if

  return
end
subroutine dtris2 ( npt, vcl, ind, ntri, til, tnbr, ierror )

!*****************************************************************************80
!
!! DTRIS2 constructs a Delaunay triangulation of 2D vertices.
!
!  Discussion:
!
!    The routine constructs the Delaunay triangulation of a set of 2D vertices
!    using an incremental approach and diagonal edge swaps.  Vertices are
!    first sorted in lexicographically increasing (X,Y) order, and
!    then are inserted one at a time from outside the convex hull.
!
!    Note that DTRIS2 or RTRIS2, the fundamental routine for constructing
!    the Delaunay triangulation, alters the input coordinate data by
!    sorting it.  This has caused me so many problems that I finally 
!    wrote a modified version of DTRIS2/RTRIS2 that undoes the sorting
!    before return.  In all other programs that use DTRIS2/RTRIS2, I
!    use the modified version, but I have left the original here in this
!    package.
!
!    On abnormal return, IERROR is set to 8, 224, or 225.
!
!  Modified:
!
!    22 July 2003
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer NPT, the number of vertices.
!
!    Input, real ( kind = 8 ) VCL(2,NPT), the coordinates of the vertices.
!
!    Input/output, integer IND(NPT), the indices in VCL of the vertices 
!    to be triangulated.  On output, IND has been permuted by the sort.
!
!    Output, integer NTRI, the number of triangles in the triangulation; 
!    NTRI is equal to 2*NPT - NB - 2, where NB is the number of boundary 
!    vertices.
!
!    Output, integer TIL(3,NTRI), the nodes that make up each triangle.
!    The elements are indices of VCL.  The vertices of the triangles are 
!    in counter clockwise order.
!
!    Output, integer TNBR(3,NTRI), the triangle neighbor list.
!    Positive elements are indices of TIL; negative elements are used for links
!    of a counter clockwise linked list of boundary edges; LINK = -(3*I + J-1)
!    where I, J = triangle, edge index; TNBR(J,I) refers to
!    the neighbor along edge from vertex J to J+1 (mod 3).
!
!    Output, integer IERROR, an error flag, nonzero if an error occurred.
!
  integer npt

  real ( kind = 8 ) cmax
  integer e
  integer i
  integer ierror
  integer ind(npt)
  integer j
  integer k
  integer l
  integer ledg
  integer lr
  integer lrline
  integer ltri
  integer m
  integer m1
  integer m2
  integer, parameter :: msglvl = 0
  integer n
  integer ntri
  integer redg
  integer rtri
  integer stack(npt)
  integer t
  integer til(3,npt*2)
  integer tnbr(3,npt*2)
  real ( kind = 8 ) tol
  integer top
  real ( kind = 8 ) vcl(2,npt)

  maxst = npt

  ierror = 0
  tol = 100.0D+00 * epsilon ( tol )

  ierror = 0
!
!  Sort the vertices.
!
  call dhpsrt ( 2, npt, 2, vcl, ind )
!
!  Ensure that no two consecutive points are too close.
!
  m1 = ind(1)

  do i = 2, npt

    m = m1
    m1 = ind(i)

    k = 0
    do j = 1, 2

      cmax = max ( abs ( vcl(j,m) ), abs ( vcl(j,m1) ) )

      if ( tol * cmax < abs ( vcl(j,m) - vcl(j,m1) ) .and. tol < cmax ) then
        k = j
        exit
      end if

    end do

    if ( k == 0 ) then
      ierror = 224
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DTRIS2 - Fatal error!'
      write ( *, '(a)' ) '  Consecutive points are too close.'
      return
    end if

  end do
!
!  Take the first two points, M1 and M2, and find a suitable non-collinear
!  third, M.  All points between M2 and M are very close to collinear
!  with M1 and M2.
!
  m1 = ind(1)
  m2 = ind(2)
  j = 3

  do

    if ( npt < j ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DTRIS2 - Fatal error!'
      write ( *, '(a)' ) '  Could not find a noncollinear third point.'
      ierror = 225
      return
    end if

    m = ind(j)
    lr = lrline ( vcl(1,m), vcl(2,m), vcl(1,m1), vcl(2,m1), vcl(1,m2), &
      vcl(2,m2), 0.0D+00 )

    if ( lr /= 0 ) then
      exit
    end if

    j = j + 1
  
  end do

  ntri = j - 2
!
!  Depending on the orientation of M1, M2, and M, set up the initial
!  triangle data.
!
  if ( lr == -1 ) then

    til(1,1) = m1
    til(2,1) = m2
    til(3,1) = m
    tnbr(3,1) = -3

    do i = 2, ntri

      m1 = m2
      m2 = ind(i+1)
      til(1,i) = m1
      til(2,i) = m2
      til(3,i) = m
      tnbr(1,i-1) = -3 * i
      tnbr(2,i-1) = i
      tnbr(3,i) = i - 1

    end do

    tnbr(1,ntri) = -3 * ntri - 1
    tnbr(2,ntri) = -5
    ledg = 2
    ltri = ntri

  else

    til(1,1) = m2
    til(2,1) = m1
    til(3,1) = m
    tnbr(1,1) = -4

    do i = 2, ntri
      m1 = m2
      m2 = ind(i+1)
      til(1,i) = m2
      til(2,i) = m1
      til(3,i) = m
      tnbr(3,i-1) = i
      tnbr(1,i) = -3 * i - 3
      tnbr(2,i) = i - 1
    end do

    tnbr(3,ntri) = -3 * ntri
    tnbr(2,1) = -3 * ntri - 2
    ledg = 2
    ltri = 1

  end if

  if ( msglvl == 4 ) then

    m2 = ind(1)
    write ( *, '(i7,4f15.7)' ) 1,vcl(1,m2),vcl(2,m2),vcl(1,m),vcl(2,m)

    do i = 2, j-1
      m1 = m2
      m2 = ind(i)
      write ( *, '(i7,4f15.7)' ) 1,vcl(1,m1),vcl(2,m1),vcl(1,m2),vcl(2,m2)
      write ( *, '(i7,4f15.7)' ) 1,vcl(1,m2),vcl(2,m2),vcl(1,m),vcl(2,m)
    end do

  end if
!
!  Insert vertices one at a time from outside the convex hull, determine
!  the visible boundary edges, and apply diagonal edge swaps until
!  the Delaunay triangulation of the vertices (so far) is obtained.
!
  top = 0

  do i = j+1, npt

    if ( msglvl == 4 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(i7)' ) i
    end if

    m = ind(i)
    m1 = til(ledg,ltri)

    if ( ledg <= 2 ) then
      m2 = til(ledg+1,ltri)
    else
      m2 = til(1,ltri)
    end if

    lr = lrline ( vcl(1,m), vcl(2,m), vcl(1,m1), vcl(2,m1), vcl(1,m2), &
      vcl(2,m2), 0.0D+00 )

    if ( 0 < lr ) then
      rtri = ltri
      redg = ledg
      ltri = 0
    else
      l = -tnbr(ledg,ltri)
      rtri = l / 3
      redg = mod ( l, 3 ) + 1
    end if

    call vbedg ( vcl(1,m), vcl(2,m), vcl, til, tnbr, ltri, ledg, rtri, redg )

    n = ntri + 1
    l = -tnbr(ledg,ltri)

    do

      t = l / 3
      e = mod ( l, 3 ) + 1
      l = -tnbr(e,t)
      m2 = til(e,t)

      if ( e <= 2 ) then
        m1 = til(e+1,t)
      else
        m1 = til(1,t)
      end if

      ntri = ntri + 1
      tnbr(e,t) = ntri
      til(1,ntri) = m1
      til(2,ntri) = m2
      til(3,ntri) = m
      tnbr(1,ntri) = t
      tnbr(2,ntri) = ntri - 1
      tnbr(3,ntri) = ntri + 1
      top = top + 1

      if ( maxst < top ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'DTRIS2 - Fatal error!'
        write ( *, '(a)' ) '  Ran out of stack space.'
        write ( *, '(a,i6)' ) '  MAXST = ', maxst
        write ( *, '(a,i6)' ) '  TOP =   ', top
        ierror = 8
        return
      end if

      stack(top) = ntri

      if ( msglvl == 4 ) then
        write ( *, '(i7,4f15.7)' ) 1,vcl(1,m),vcl(2,m),vcl(1,m2),vcl(2,m2)
      end if

      if ( t == rtri .and. e == redg ) then
        exit
      end if

    end do

    if ( msglvl == 4 ) then
      write ( *, '(i7,4f15.7)' ) 1,vcl(1,m),vcl(2,m),vcl(1,m1),vcl(2,m1)
    end if

    tnbr(ledg,ltri) = -3 * n - 1
    tnbr(2,n) = -3 * ntri - 2
    tnbr(3,ntri) = -l
    ltri = n
    ledg = 2

    call swapec ( m, top, maxst, ltri, ledg, vcl, til, tnbr, stack, ierror )

    if ( ierror /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DTRIS2 - Fatal error!'
      write ( *, '(a,i6)' ) '  SWAPEC returned IERROR = ', ierror
      return
    end if

  end do

  if ( msglvl == 4 ) then
    write ( *, '(i7)' ) npt + 1
  end if

  return
end
subroutine dtriw2 ( npt, maxst, vcl, ind, ntri, til, tnbr, stack, ierror )

!*****************************************************************************80
!
!! DTRIW2 constructs an incremental Delaunay triangulation in 2D.
!
!  Purpose: 
!
!    Construct the Delaunay triangulation of 2D vertices using
!    incremental approach and diagonal edge swaps. Vertices are
!    inserted one at a time in order given by IND array. The initial
!    triangles created due to a new vertex are obtained by a walk
!    through the triangulation until location of vertex is known.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer NPT, the number of 2D points (vertices).
!
!    Input, integer MAXST, the maximum size available for STACK array; should 
!    be about NPT to be safe, but MAX(10,2*LOG2(NPT)) usually enough.
!
!    Input, real ( kind = 8 ) VCL(1:2,1:*), the coordinates of 2D vertices.
!
!    Input, integer IND(1:NPT), indices in VCL of vertices to be triangulated;
!    vertices are inserted in order given by this array.
!
!    Output, integer NTRI, the number of triangles in triangulation; equal to
!    2*NPT - NB - 2 where NB = number of boundary vertices.
!
!    Output, integer TIL(1:3,1:NTRI), the triangle incidence list; elements 
!    are indices of VCL; vertices of triangles are in counter clockwise order.
!
!    Output, integer TNBR(1:3,1:NTRI), the triangle neighbor list; positive 
!    elements are indices of TIL; negative elements are used for links
!    of counter clockwise linked list of boundary edges; LINK = -(3*I + J-1)
!    where I, J = triangle, edge index; TNBR(J,I) refers to
!    the neighbor along edge from vertex J to J+1 (mod 3).
!
!    Workspace, integer STACK(1:MAXST), used for stack of triangles for which
!    circumcircle test must be made.
!
!    Output, integer IERROR, error flag. For abnormal return,
!    IERROR is set to 8, 224, 225, or 226.
!
  integer maxst
  integer npt

  integer bedg
  integer btri
  real ( kind = 8 ) cmax
  integer e
  integer em1
  integer ep1
  integer ntri
  integer i
  integer i3
  integer ierror
  integer ind(npt)
  integer j
  integer l
  integer ledg
  integer lr
  integer lrline
  integer ltri
  integer m
  integer m1
  integer m2
  integer m3
  integer, parameter :: msglvl = 0
  integer n
  integer redg
  integer rtri
  integer stack(maxst)
  integer t
  integer til(3,npt*2)
  integer tnbr(3,npt*2)
  integer top
  real ( kind = 8 ) tol
  real ( kind = 8 ) vcl(2,*)

  ierror = 0
  tol = 100.0D+00 * epsilon ( tol )
!
!  Determine the initial triangle.
!
  m1 = ind(1)
  m2 = ind(2)

  do j = 1, 2
    cmax = max ( abs ( vcl(j,m1) ), abs ( vcl(j,m2) ) )
    if ( tol * cmax < abs ( vcl(j,m1) - vcl(j,m2) ) .and. tol < cmax ) then
      go to 20
    end if
  end do

  ierror = 224
  return

20 continue

  i3 = 3

30 continue

  if ( npt < i3 ) then
    ierror = 225
    return
  end if

  m = ind(i3)
  lr = lrline ( vcl(1,m), vcl(2,m), vcl(1,m1), vcl(2,m1), vcl(1,m2), &
    vcl(2,m2), 0.0D+00 )

  if ( lr == 0 ) then
    i3 = i3 + 1
    go to 30
  end if

  if ( i3 /= 3 ) then
    ind(i3) = ind(3)
    ind(3) = m
  end if

  ntri = 1

  if ( lr == -1 ) then
    til(1,1) = m1
    til(2,1) = m2
  else
    til(1,1) = m2
    til(2,1) = m1
  end if

  til(3,1) = m
  tnbr(1,1) = -4
  tnbr(2,1) = -5
  tnbr(3,1) = -3

  if ( msglvl == 4 ) then
    write ( *,600) 1,vcl(1,m1),vcl(2,m1),vcl(1,m2),vcl(2,m2)
    write ( *,600) 1,vcl(1,m2),vcl(2,m2),vcl(1,m),vcl(2,m)
    write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m1),vcl(2,m1)
  end if
!
!  Insert vertices one at a time from anywhere.
!  Walk through the triangulation to determine the location of the new vertex.
!  Apply diagonal edge swaps until Delaunay triangulation of vertices
!  (so far) is obtained.
!
  top = 0

  do i = 4, npt

    if ( msglvl == 4 ) then
      write ( *,600) i
    end if

    m = ind(i)
    rtri = ntri

    call walkt2 ( vcl(1,m), vcl(2,m), ntri, vcl, til, tnbr, rtri, redg, ierror )

    if ( redg == 0 ) then

      m1 = til(1,rtri)
      m2 = til(2,rtri)
      m3 = til(3,rtri)
      til(3,rtri) = m

      if ( 0 < tnbr(1,rtri) ) then
        top = 1
        stack(top) = rtri
      end if

      ntri = ntri + 1
      til(1,ntri) = m2
      til(2,ntri) = m3
      til(3,ntri) = m
      n = tnbr(2,rtri)
      tnbr(1,ntri) = n

      if ( 0 < n ) then
        if ( tnbr(1,n) == rtri ) then
          tnbr(1,n) = ntri
        else if ( tnbr(2,n) == rtri ) then
          tnbr(2,n) = ntri
        else
          tnbr(3,n) = ntri
        end if
        top = top + 1
        stack(top) = ntri
      end if

      ntri = ntri + 1
      til(1,ntri) = m3
      til(2,ntri) = m1
      til(3,ntri) = m
      n = tnbr(3,rtri)
      tnbr(1,ntri) = n

      if ( 0 < n ) then
        if ( tnbr(1,n) == rtri ) then
          tnbr(1,n) = ntri
        else if ( tnbr(2,n) == rtri ) then
          tnbr(2,n) = ntri
        else
          tnbr(3,n) = ntri
        end if
        top = top + 1
        stack(top) = ntri
      end if

      tnbr(2,rtri) = ntri - 1
      tnbr(3,rtri) = ntri
      tnbr(2,ntri-1) = ntri
      tnbr(3,ntri-1) = rtri
      tnbr(2,ntri) = rtri
      tnbr(3,ntri) = ntri - 1

      if ( tnbr(1,ntri-1) <= 0 ) then

        t = rtri
        e = 1

        do

          if ( tnbr(e,t) <= 0 ) then
            exit
          end if

          t = tnbr(e,t)

          if ( til(1,t) == m2 ) then
            e = 3
          else if ( til(2,t) == m2 ) then
            e = 1
          else
            e = 2
          end if

        end do

        tnbr(e,t) = -3 * ntri + 3

      end if

      if ( tnbr(1,ntri) <= 0 ) then

        t = ntri - 1
        e = 1

        do

          if ( tnbr(e,t) <= 0 ) then
            exit
          end if

          t = tnbr(e,t)
          if ( til(1,t) == m3 ) then
            e = 3
          else if ( til(2,t) == m3 ) then
            e = 1
          else
            e = 2
          end if

        end do

        tnbr(e,t) = -3 * ntri

      end if

      if ( msglvl == 4 ) then
        write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m1),vcl(2,m1)
        write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m2),vcl(2,m2)
        write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m3),vcl(2,m3)
      end if

   else if ( redg < 0 ) then

      redg = -redg
      ltri = 0
      call vbedg ( vcl(1,m), vcl(2,m), vcl, til, tnbr, ltri, ledg, rtri, redg )
      n = ntri + 1
      l = -tnbr(ledg,ltri)

60    continue

         t = l / 3
         e = mod ( l, 3 ) + 1
         l = -tnbr(e,t)
         m2 = til(e,t)

         if ( e <= 2 ) then
           m1 = til(e+1,t)
         else
           m1 = til(1,t)
         end if

         ntri = ntri + 1
         tnbr(e,t) = ntri
         til(1,ntri) = m1
         til(2,ntri) = m2
         til(3,ntri) = m
         tnbr(1,ntri) = t
         tnbr(2,ntri) = ntri - 1
         tnbr(3,ntri) = ntri + 1
         top = top + 1

         if ( maxst < top ) then
           ierror = 8
           go to 100
         end if
 
        stack(top) = ntri

         if ( msglvl == 4 ) then
           write (*,600) 1,vcl(1,m),vcl(2,m),vcl(1,m2),vcl(2,m2)
         end if

      if ( t /= rtri .or. e /= redg ) then
        go to 60
      end if

      if ( msglvl == 4 ) then
        write (*,600) 1,vcl(1,m),vcl(2,m),vcl(1,m1),vcl(2,m1)
      end if

      tnbr(ledg,ltri) = -3*n - 1
      tnbr(2,n) = -3*ntri - 2
      tnbr(3,ntri) = -l

   else if ( redg <= 3 ) then

      m1 = til(redg,rtri)

      if ( redg == 1 ) then
        e = 2
        ep1 = 3
      else if ( redg == 2 ) then
        e = 3
        ep1 = 1
      else
        e = 1
        ep1 = 2
      end if

      m2 = til(e,rtri)
      til(e,rtri) = m
      m3 = til(ep1,rtri)

      if ( 0 < tnbr(ep1,rtri) ) then
        top = 1
        stack(top) = rtri
      end if

      ntri = ntri + 1
      til(1,ntri) = m
      til(2,ntri) = m2
      til(3,ntri) = m3
      n = tnbr(e,rtri)
      tnbr(2,ntri) = n
      tnbr(3,ntri) = rtri
      tnbr(e,rtri) = ntri

      if ( 0 < n ) then
        if ( tnbr(1,n) == rtri ) then
          tnbr(1,n) = ntri
        else if ( tnbr(2,n) == rtri ) then
          tnbr(2,n) = ntri
        else
          tnbr(3,n) = ntri
        end if
        top = top + 1
        stack(top) = ntri
      end if

      if ( msglvl == 4 ) then
        write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m3),vcl(2,m3)
      end if

      ltri = tnbr(redg,rtri)

      if ( ltri <= 0 ) then
        tnbr(1,ntri) = ltri
        tnbr(redg,rtri) = -3*ntri
        if ( tnbr(2,ntri) <= 0 ) then
          tnbr(1,ntri) = -3*ntri - 1
        end if

      else

        tnbr(1,ntri) = ntri + 1
        tnbr(redg,rtri) = ltri

        if ( til(1,ltri) == m2 ) then
          ledg = 1
          em1 = 2
          e = 3
        else if ( til(2,ltri) == m2 ) then
          ledg = 2
          em1 = 3
          e = 1
        else
          ledg = 3
          em1 = 1
          e = 2
        end if

        til(ledg,ltri) = m
        m3 = til(e,ltri)
        if ( 0 < tnbr(em1,ltri) ) then
          top = top + 1
          stack(top) = ltri
        end if
        ntri = ntri + 1
        til(1,ntri) = m2
        til(2,ntri) = m
        til(3,ntri) = m3
        tnbr(1,ntri) = ntri - 1
        tnbr(2,ntri) = ltri
        n = tnbr(e,ltri)
        tnbr(3,ntri) = n
        tnbr(e,ltri) = ntri

        if ( 0 < n ) then
          if ( tnbr(1,n) == ltri ) then
            tnbr(1,n) = ntri
          else if ( tnbr(2,n) == ltri ) then
            tnbr(2,n) = ntri
          else
            tnbr(3,n) = ntri
          end if
          top = top + 1
          stack(top) = ntri
        end if

        if ( msglvl == 4 ) then
          write ( *,600) 1,vcl(1,m),vcl(2,m),vcl(1,m3),vcl(2,m3)
        end if

        if ( tnbr(2,ntri-1) <= 0 ) then

          t = ntri
          e = 3

          do

            if ( tnbr(e,t) <= 0 ) then
              exit
            end if

            t = tnbr(e,t)
            if ( til(1,t) == m2 ) then
              e = 3
            else if ( til(2,t) == m2 ) then
              e = 1
            else
              e = 2
            end if

          end do

          tnbr(e,t) = -3 * ntri + 2

        end if

         if ( tnbr(3,ntri) <= 0 ) then

          t = ltri

          if ( ledg <= 2 ) then
             e = ledg + 1
          else
             e = 1
          end if

          do

            if ( tnbr(e,t) <= 0 ) then
              exit
            end if

            t = tnbr(e,t)
            if ( til(1,t) == m3 ) then
              e = 3
            else if ( til(2,t) == m3 ) then
              e = 1
            else
              e = 2
            end if

          end do

          tnbr(e,t) = -3 * ntri - 2

         end if

      end if

    else
      ierror = 224
      go to 100
    end if

    btri = 0
    bedg = 0

    call swapec ( m, top, maxst, btri, bedg, vcl, til, tnbr, stack, ierror )

    if ( ierror /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DTRIW2 - Fatal error!'
      write ( *, '(a,i6)' ) '  SWAPEC returned error flag IERROR = ', ierror
      exit
    end if

  end do

100 continue

  if ( i3 /= 3 ) then
    t = ind(i3)
    ind(i3) = ind(3)
    ind(3) = t
  end if

  if ( msglvl == 4 ) then
    write ( *,600) npt+1
  end if

600 format (1x,i7,4f15.7)

  return
end
subroutine edght ( a, b, v, n, htsiz, maxedg, hdfree, last, ht, edge, w, &
  ierror )

!*****************************************************************************80
!
!! EDGHT searches a hash table for a record in EDGE containing key (A,B).
!
!  Purpose: 
!
!    Search in hash table HT for record in EDGE containing
!    key (A,B).
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer A, B, vertex indices, greater than 0, of edge (also 
!    key of hash table).
!
!    Input, integer V, value associated with edge.
!
!    Input, integer N, upper bound on A, B.
!
!    Input, integer HTSIZ, the size of hash table HT.
!
!    Input, integer MAXEDG, the maximum size available for EDGE array.
!
!    Input/output, integer HDFREE, head pointer to linked list of free entries 
!    of EDGE array due to deletions.  Before first call to this routine, HDFREE
!    should be set to 0.
!
!    Input/output, integer LAST, index of last entry used in EDGE array.  
!    Before first call to this routine, LAST should be set to 0.
!
!    Input/output, integer HT(0:HTSIZ-1), hash table of head pointers (direct 
!    chaining with ordered lists is used).  Before first call to this routine, 
!    entries of HT should be set to 0.  If key with A,B is found then this 
!    record is deleted from hash table, else record is inserted in hash table.
!
!    Input/output, integer EDGE(1:4,1:MAXEDG), entries of hash table records;
!      EDGE(1,I) = MIN(A,B); EDGE(2,I) = MAX(A,B);
!      EDGE(3,I) = V; EDGE(4,I) = link
!    If key with A,B is found then this record is deleted
!    from hash table, else record is inserted in hash table.
!
!    Output, integer W, EDGE(3,INDEX), where INDEX is index of record, if found;
!    else 0.
!
!    Output, integer IERROR, error flag.  For abnormal return,
!    IERROR is set to 1
!
  integer htsiz
  integer maxedg

  integer a
  integer aa
  integer b
  integer bb
  integer bptr
  integer edge(4,maxedg)
  integer hdfree
  integer ht(0:htsiz-1)
  integer ierror
  integer k
  integer last
  integer n
  integer newp
  integer ptr
  integer v
  integer w

  ierror = 0

  if ( a < b ) then
    aa = a
    bb = b
  else
    aa = b
    bb = a
  end if

  k = mod ( aa*n + bb, htsiz )
  bptr = -1
  ptr = ht(k)

10 continue

  if ( ptr /= 0 ) then

    if ( aa < edge(1,ptr) ) then

      go to 20

    else if ( edge(1,ptr) == aa ) then

      if ( bb < edge(2,ptr) ) then

        go to 20

      else if ( edge(2,ptr) == bb ) then

        if ( bptr == -1 ) then
          ht(k) = edge(4,ptr)
        else
          edge(4,bptr) = edge(4,ptr)
        end if

        edge(4,ptr) = hdfree
        hdfree = ptr
        w = edge(3,ptr)
        return

      end if

    end if

    bptr = ptr
    ptr = edge(4,ptr)
    go to 10

  end if

20 continue

  if ( 0 < hdfree ) then

    newp = hdfree
    hdfree = edge(4,hdfree)

  else

    last = last + 1
    newp = last

    if ( maxedg < last ) then
      ierror = 1
      return
    end if

  end if

  if ( bptr == -1 ) then
    ht(k) = newp
  else
    edge(4,bptr) = newp
  end if

  edge(1,newp) = aa
  edge(2,newp) = bb
  edge(3,newp) = v
  edge(4,newp) = ptr
  w = 0

  return
end
subroutine eqdis2 ( hflag, umdf, kappa, angspc, angtol, dmin, nmin, ntrid,  &
  nvc, npolg, nvert, maxvc, maxhv, maxpv, maxiw, maxwk, vcl, regnum, hvl,  &
  pvl, iang, area, psi, h, iwk, wk, ierror )

!*****************************************************************************80
!
!! EQDIS2 further subdivides convex polygons for mesh equidistribution.
!
!  Purpose: 
!
!    Further subdivide convex polygons so that an approximately
!    equidistributing triangular mesh can be constructed with
!    respect to a heuristic or a user-supplied mesh distribution
!    function (MDF), and determine triangle size for each polygon of
!    decomposition.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, logical HFLAG, is .TRUE. if heuristic mdf, .FALSE. if
!    user-supplied mdf.
!
!    Input, external UMDF, a user-supplied mesh distribution function 
!    of the form:
!
!      function umdf ( x, y )
!      double precision umdf
!      double precision x
!      double precision y
!
!    Input, real ( kind = 8 ) KAPPA, the mesh smoothness parameter in 
!    the interval [0.0,1.0].
!
!    Input, real ( kind = 8 ) ANGSPC, the angle spacing parameter in radians
!    used to determine extra points as possible endpoints of separators.
!
!    Input, real ( kind = 8 ) ANGTOL, the angle tolerance parameter in radians
!    used in accepting separators.
!
!    Input, real ( kind = 8 ) DMIN, a parameter used to determine if variation 
!    of mdf in polygon is 'sufficiently high'.
!
!    Input, integer NMIN, a parameter used to determine if 'sufficiently large'
!    number of triangles in polygon.
!
!    Input, integer NTRID, the desired number of triangles in mesh.
!
!    Input/output, integer NVC, the number of vertex coordinates or 
!    positions used in VCL array.
!
!    Input/output, integer NPOLG, the number of polygonal subregions or 
!    positions used in HVL array.
!
!    Input/output, integer NVERT, the number of polygon vertices or positions 
!    used in PVL array.
!
!    Input, integer MAXVC, the maximum size available for VCL array, should 
!    be greater than or equal to the number of vertex coordinates required
!    for decomposition (approximately NVC + 2*NS where NS is expected number 
!    of new separators).
!
!    Input, integer MAXHV, the maximum size available for HVL, REGNUM, AREA, 
!    PSI, H arrays; should be greater than or equal to the number of polygons
!    required for decomposition (approximately NPOLG + NS).
!
!    Input, integer MAXPV, the maximum size available for PVL, IANG arrays; 
!    should be greater than or equal to the number of polygon vertices 
!    required for decomposition (approximately NVERT + 5*NS).
!
!    Input, integer MAXIW, the maximum size available for IWK array; should 
!    be greater than or equal to 
!      MAX(2*NP, NVERT + NPOLG + 3*NVRT + INT(2*PI/ANGSPC))
!    where NVRT is maximum number of vertices in a convex
!    polygon of the (input) decomposition, NP is expected
!    value of NPOLG on output.
!
!    Input, integer MAXWK, the maximum size available for WK array; should
!    be greater than or equal to 
!      NVC + NVERT + 2*NPOLG + 3*(NVRT + INT(2*PI/ANGSPC)).
!
!    Input/output, real ( kind = 8 ) VCL(1:2,1:NVC), the vertex coordinate list.
!
!    Input/output, integer REGNUM(1:NPOLG), the region numbers.
!
!    Input/output, real ( kind = 8 ) HVL(1:NPOLG), head vertex list.
!
!    Input/output, integer PVL(1:4,1:NVERT), real ( kind = 8 ) IANG(1:NVERT), 
!    the polygon vertex list and interior angles; see routine DSPGDC for more 
!    details.  Note that the data structures should be as output from routine
!    CVDEC2.
!
!    Output, real ( kind = 8 ) AREA(1:NPOLG), the area of convex polygons 
!    in the decomposition.
!
!    Output, real ( kind = 8 ) PSI(1:NPOLG), the smoothed mean mdf 
!    values in the convex polygons.
!
!    Output, real ( kind = 8 ) H(1:NPOLG), the triangle size for convex
!    polygons.
!
!    Workspace, integer IWK(1:MAXIW).
!
!    Workspace, real ( kind = 8 ) WK(1:MAXWK).
!
!    Output, integer IERROR, error flag.  For abnormal return,
!    IERROR is set to 3, 4, 5, 6, 7, 200, 201, or 222.
!
  integer maxhv
  integer maxiw
  integer maxpv
  integer maxvc
  integer maxwk

  real ( kind = 8 ) angspc
  real ( kind = 8 ) angtol
  real ( kind = 8 ) area(maxhv)
  real ( kind = 8 ) dmin
  integer edgval
  real ( kind = 8 ) h(maxhv)
  logical hflag
  integer hvl(maxhv)
  real ( kind = 8 ) iang(maxpv)
  integer ierror
  integer ivrt
  integer iwk(maxiw)
  real ( kind = 8 ) kappa
  integer m
  integer n
  integer nmin
  integer npolg
  integer ntrid
  integer nvc
  integer nvert
  real ( kind = 8 ) psi(maxhv)
  integer pvl(4,maxpv)
  integer regnum(maxhv)
  real ( kind = 8 ) umdf
  external umdf
  real ( kind = 8 ) vcl(2,maxvc)
  integer vrtval
  integer widsq
  real ( kind = 8 ) wk(maxwk)
  integer xivrt

  ierror = 0
  ivrt = 1
  xivrt = ivrt + nvert
  m = xivrt + npolg

  if ( maxiw < m ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EQDIS2 - Fatal error!'
    write ( *, '(a)' ) '  MAXIW < M'
    ierror = 6
    return
  end if

  widsq = 1

  if ( hflag ) then

    edgval = widsq + npolg
    vrtval = edgval + nvert
    n = npolg + nvert + nvc

    if ( maxwk < n ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'EQDIS2 - Fatal error!'
      write ( *, '(a)' ) '  MAXWK < N'
      ierror = 7
      return
    end if

  else

    edgval = 1
    vrtval = 1
    n = 0

  end if

  call dsmdf2 ( hflag, nvc, npolg, maxwk-n, vcl, hvl, pvl, iang, iwk(ivrt), &
    iwk(xivrt), wk(widsq), wk(edgval), wk(vrtval), area, wk(n+1), ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EQDIS2 - Fatal error!'
    write ( *, '(a)' ) '  DSMDF2 returns error condition.'
    return
  end if

  call mfdec2 ( hflag, umdf, kappa, angspc, angtol, dmin, nmin, ntrid, nvc, &
    npolg, nvert, maxvc, maxhv, maxpv, maxiw-m, maxwk-n, vcl, regnum, hvl, &
    pvl, iang, iwk(ivrt), iwk(xivrt), wk(widsq), wk(edgval), wk(vrtval), &
    area, psi, iwk(m+1), wk(n+1), ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EQDIS2 - Fatal error!'
    write ( *, '(a)' ) '  MFDEC2 returns error condition.'
    return
  end if

  if ( maxiw < 2 * npolg ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'EQDIS2 - Fatal error!'
    write ( *, '(a)' ) '  MAXIW < 2 * NPOLG.'
    ierror = 6
    return
  end if

  call trisiz ( ntrid, npolg, hvl, pvl, area, psi, h, iwk, iwk(npolg+1) )

  return
end
subroutine fndsep ( angac1, xr, yr, nvrt, xc, yc, ivis, theta, nv, iv,  &
  vcl, pvl, iang, angsep, i1, i2, wkang )

!*****************************************************************************80
!
!! FNDSEP finds separators to resolve a reflex vertex.
!
!  Purpose: 
!
!    Find 1 or 2 separators which can resolve a reflex vertex
!    (XR,YR) using a max-min angle criterion from list of vertices
!    in increasing polar angle with respect to the reflex vertex.
!
!    Preference is given to 1 separator.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ANGAC1, the angle tolerance parameter used 
!    for preference in accepting one separator.
!
!    Input, real ( kind = 8 ) XR, YR, the coordinates of reflex vertex.
!
!    Input, integer NVRT, (number of vertices) - 1.
!
!    Input, real ( kind = 8 ) XC(0:NVRT), YC(0:NVRT), the vertex coordinates 
!    of possible endpoints of a separator.
!
!    Input, integer IVIS(0:NVRT), contains information about the vertices of
!    XC, YC arrays with respect to the polygon vertex list; if
!    0 < IVIS(I) then vertex (XC(I),YC(I)) has index IVIS(I)
!    in PVL; if IVIS(I) < 0 then vertex (XC(I),YC(I)) is on
!    the edge joining vertices with indices -IVIS(I) and
!    SUCC(-IVIS(I)) in PVL.
!
!    Input, real ( kind = 8 ) THETA(0:NVRT), the polar angles of vertices 
!    in increasing order; THETA(NVRT) is the interior angle of reflex vertex;
!    THETA(I), 0 <= I, is the polar angle of (XC(I),YC(I))
!    with respect to reflex vertex.
!
!    Input, integer NV, (number of vertices to be considered as endpoint of a
!    separator) - 1.
!
!    Input, integer IV(0:NV), the indices of vertices in XC, YC arrays to be
!    considered as endpoint of a separator; angle between
!    consecutive vertices is assumed to be < 180 degrees.
!
!    Input, real ( kind = 8 ) VCL(1:2,1:*), the vertex coordinate list.
!
!    Input, integer PVL(1:4,1:*), real ( kind = 8 ) IANG(1:*), the polygon 
!    vertex list, interior angles.
!
!    Output, real ( kind = 8 ) ANGSEP, the minimum of the 4 or 7 angles at the
!    boundary resulting from 1 or 2 separators, respectively.
!
!    Output, integer I1, I2, the indices of endpoints of separators in XC, 
!    YC arrays; I2 = -1 if there is only one separator, else I1 < I2.
!
!    Workspace, real ( kind = 8 ) WKANG(0:NV).
!
  real ( kind = 8 ) ang
  real ( kind = 8 ) angac1
  real ( kind = 8 ) angsep
  real ( kind = 8 ) angsp2
  integer i
  integer i1
  integer i2
  real ( kind = 8 ) iang(*)
  integer ii
  integer k
  integer l
  integer m
  integer nl
  integer nr
  integer nv
  integer nvrt
  integer iv(0:nv)
  integer ivis(0:nvrt)
  real ( kind = 8 ) minang
  integer p
  real ( kind = 8 ) phi
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  integer pvl(4,*)
  integer q
  integer r
  real ( kind = 8 ) theta(0:nvrt)
  real ( kind = 8 ) tol
  real ( kind = 8 ) vcl(2,*)
  real ( kind = 8 ) wkang(0:nv)
  real ( kind = 8 ) xc(0:nvrt)
  real ( kind = 8 ) xr
  real ( kind = 8 ) yc(0:nvrt)
  real ( kind = 8 ) yr

  tol = 100.0D+00 * epsilon ( tol )
!
!  Determine the vertices in the inner cone - indices P to Q.
!
  i = 0
  p = -1
  phi = theta(nvrt) - pi + tol

  do while ( p < 0 )

    if ( phi <= theta(iv(i)) ) then
      p = i
    else
      i = i + 1
    end if

  end do

  i = nv
  q = -1
  phi = pi - tol

  do while ( q < 0 )

    if ( theta(iv(i)) <= phi ) then
      q = i
    else
      i = i - 1
    end if

  end do
!
!  Use the max-min angle criterion to find the best separator
!  in inner cone.
!
  angsep = 0.0

  do i = p, q

    k = iv(i)
    ang = minang ( xr, yr, xc(k), yc(k), ivis(k), theta(k), theta(nvrt), &
      vcl, pvl, iang )

    if ( angsep < ang ) then
      angsep = ang
      ii = iv(i)
    end if

  end do

  angsp2 = angsep
  if ( angac1 <= angsep ) then
    go to 110
  end if
!
!  If the best separator in inner cone is not 'good' enough,
!  use max-min angle criterion to try to find a better pair
!  of separators from the right and left cones.
!
  nr = 0
  nl = 0

  do r = 0, p-1

    wkang(r) = 0.0D+00

    if ( angsep < theta(iv(r)) ) then

      k = iv(r)

      ang = minang ( xr, yr, xc(k), yc(k), ivis(k), theta(k), theta(nvrt), &
        vcl, pvl, iang )

      if ( angsep < ang ) then
        nr = nr + 1
        wkang(r) = ang
      end if

    end if

  end do

  if ( nr == 0 ) then
    go to 110
  end if

  phi = theta(nvrt) - angsep

  do l = q+1, nv

    wkang(l) = 0.0D+00

    if ( theta(iv(l)) < phi ) then

      k = iv(l)
      ang = minang ( xr, yr, xc(k), yc(k), ivis(k), theta(k), theta(nvrt), &
        vcl, pvl, iang )

      if ( angsep < ang ) then
        nl = nl + 1
        wkang(l) = ang
      end if

    end if

  end do

  if ( nl == 0 ) then
    go to 110
  end if
!
!  Check all possible pairs for the best pair of separators
!  in the right and left cones.
!
  m = nv

  do r = p-1, 0, -1

     if ( q < m .and. angsp2 < wkang(r) ) then

      phi = theta(iv(r))

80    continue

      if ( q < m .and. &
           ( wkang(m) <= angsp2 .or. &
             pi - tol < theta(iv(m)) - phi ) ) then
        m = m - 1
        go to 80
      end if

      do l = q+1, m

         if ( angsp2 < wkang(l) ) then

            ang = min ( theta(iv(l)) - phi, wkang(r), wkang(l) )

            if ( angsp2 < ang ) then
             angsp2 = ang
             i1 = iv(r)
             i2 = iv(l)
            end if

         end if

       end do

     end if

  end do
!
!  Choose 1 or 2 separators based on max-min angle criterion or
!  ANGAC1 parameter.
!
  110 continue

  if ( angsp2 <= angsep ) then
    i1 = ii
    i2 = -1
  else
    angsep = angsp2
  end if

  return
end
subroutine fndtri ( iedg, mxtr, sflag, tedg, itr, ind, ierror )

!*****************************************************************************80
!
!! FNDTRI finds two triangles containing a given edge.
!
!  Purpose: 
!
!    Find two triangles containing edge with index IEDG in array TEDG.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer IEDG, the index of edge to be searched in TEDG.
!
!    Input, integer MXTR, the maximum index of triangle to be searched in TEDG.
!
!    Input, logical SFLAG, is .TRUE. if and only if the second triangle is to be 
!    searched from end of array.
!
!    Input, integer TEDG(1:3,1:MXTR), triangle edge indices; see routine CVDTRI.
!
!    Output, integer ITR(1:2), IND(1:2), indices such that IEDG =
!    TEDG(IND(1),ITR(1)) = TEDG(IND(2),ITR(2)).
!
!    Output, integer IERROR, error flag.  On abnormal return,
!    IERROR is set to 231.
!
  integer mxtr

  integer i
  integer iedg
  integer ierror
  integer ind(2)
  integer itr(2)
  integer j
  integer k
  logical sflag
  integer tedg(3,mxtr)
!
!  Search from end of array TEDG.
!
  ierror = 0
  k = 1
  j = 1
  i = mxtr

10 continue

  do

    if ( tedg(j,i) == iedg ) then
      exit
    end if

    j = j + 1

    if ( 3 < j ) then
      j = 1
      i = i - 1
      if ( i <= 0 ) then
        ierror = 231
        return
      end if
    end if

  end do

  itr(k) = i
  ind(k) = j

  if ( k == 2 ) then
    return
  end if

  k = 2

  if ( sflag ) then

    j = 1
    i = i - 1

    if ( i <= 0 ) then
      ierror = 231
      return
    end if

    go to 10

  end if
!
!  Search from beginning of array TEDG for second triangle.
!
  j = 1
  i = 1

20 continue

  if ( itr(1) <= i ) then
    ierror = 231
    return
  end if

30 continue

  if ( tedg(j,i) /= iedg ) then
    j = j + 1
    if ( 3 < j ) then
      j = 1
      i = i + 1
      go to 20
    else
      go to 30
    end if
  end if

  itr(2) = i
  ind(2) = j

  return
end
subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5, 6 and 9, which
!    are commonly reserved for console I/O).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
!  Modified:
!
!    18 September 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT, the free unit number.
!
  implicit none

  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
subroutine gtime ( time )

!*****************************************************************************80
!
!! GTIME gets the current CPU time in seconds.
!
!  Modified:
!
!    05 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real ( kind = 8 ) TIME, the current reading of the CPU clock
!    in seconds.
!
  integer clock_count
  integer clock_max
  integer clock_rate
  real ( kind = 8 ) time

  call system_clock ( clock_count, clock_rate, clock_max )

  time = real ( clock_count, kind = 8 ) / real ( clock_rate, kind = 8 )

  return
end
subroutine holvrt ( nhole, vcl, hvl, pvl, holv )

!*****************************************************************************80
!
!! HOLVRT determines top and bottom vertices of holes in polygonal regions.
!
!  Purpose: 
!
!    Determine top and bottom vertices of holes in polygonal
!    regions, and sort top vertices in decreasing (y,x) order
!    and bottom vertices in increasing (y,x) order.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer NHOLE, the number of holes in region(s).
!
!    Input, real ( kind = 8 ) VCL(1:2,1:*), the vertex coordinate list.
!
!    Input, integer HVL(1:NHOLE), the head vertex list; HVL(I) is index in 
!    PVL of head vertex of Ith hole.
!
!    Input, integer PVL(1:4,1:*), the polygon vertex list; see routine DSPGDC.
!
!    Output, integer HOLV(1:NHOLE*2), the indices in PVL of top and bottom 
!    vertices of holes; first (last) NHOLE entries are for top (bottom)
!    vertices; top (bottom) vertices are sorted in decreasing
!    (increasing) lexicographic (y,x) order of coordinates.
!
  integer nhole

  integer holv(nhole*2)
  integer hv
  integer hvl(nhole)
  integer i
  integer imax
  integer imin
  integer iv
  integer j
  integer, parameter :: loc = 1
  integer lv
  integer nhp1
  integer pvl(4,*)
  integer, parameter :: succ = 3
  real ( kind = 8 ) vcl(2,*)
  real ( kind = 8 ) x
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) y
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin
!
!  Determine top and bottom vertices of holes.
!
  do i = 1, nhole

    hv = hvl(i)
    iv = hv

    do

      lv = pvl(loc,iv)

      if ( iv == hv ) then

        imin = iv
        imax = iv
        xmin = vcl(1,lv)
        ymin = vcl(2,lv)
        xmax = xmin
        ymax = ymin

      else

        x = vcl(1,lv)
        y = vcl(2,lv)

        if ( y < ymin .or. y == ymin .and. x < xmin ) then
          imin = iv
          xmin = x
          ymin = y
        else if ( ymax < y .or. y == ymax .and. xmax < x ) then
          imax = iv
          xmax = x
          ymax = y
        end if

      end if

      iv = pvl(succ,iv)

      if ( iv == hv) then
        exit
      end if

    end do

    holv(i) = imax
    holv(i+nhole) = imin

  end do
!
!  Use linear insertion sort to sort the top vertices of the holes
!  in decreasing (y,x) order, then bottom vertices in increasing
!  (y,x) order.  It is assumed that NHOLE is small.
!
  do i = 2, nhole

    hv = holv(i)
    lv = pvl(loc,hv)
    x = vcl(1,lv)
    y = vcl(2,lv)
    j = i

30  continue

    iv = holv(j-1)
    lv = pvl(loc,iv)

    if ( vcl(2,lv) < y .or. ( y == vcl(2,lv) .and. vcl(1,lv) < x ) ) then
      holv(j) = iv
      j = j - 1
      if ( 1 < j ) then
        go to 30
      end if
    end if

    holv(j) = hv

  end do

  nhp1 = nhole + 1

  do i = nhp1+1, nhole+nhole

    hv = holv(i)
    lv = pvl(loc,hv)
    x = vcl(1,lv)
    y = vcl(2,lv)
    j = i

50  continue

    iv = holv(j-1)
    lv = pvl(loc,iv)

    if ( y < vcl(2,lv) .or. y == vcl(2,lv) .and. x < vcl(1,lv) ) then

      holv(j) = iv
      j = j - 1

      if ( nhp1 < j ) then
        go to 50
      end if

    end if

    holv(j) = hv

  end do

  return
end
function i4_modp ( i, j )

!*****************************************************************************80
!
!! I4_MODP returns the nonnegative remainder of integer division.
!
!  Formula:
!
!    If
!      NREM = I4_MODP ( I, J )
!      NMULT = ( I - NREM ) / J
!    then
!      I = J * NMULT + NREM
!    where NREM is always nonnegative.
!
!  Comments:
!
!    The MOD function computes a result with the same sign as the
!    quantity being divided.  Thus, suppose you had an angle A,
!    and you wanted to ensure that it was between 0 and 360.
!    Then mod(A,360) would do, if A was positive, but if A
!    was negative, your result would be between -360 and 0.
!
!    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
!
!  Examples:
!
!        I     J     MOD  I4_MODP    Factorization
!
!      107    50       7       7    107 =  2 *  50 + 7
!      107   -50       7       7    107 = -2 * -50 + 7
!     -107    50      -7      43   -107 = -3 *  50 + 43
!     -107   -50      -7      43   -107 =  3 * -50 + 43
!
!  Modified:
!
!    02 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, the number to be divided.
!
!    Input, integer J, the number that divides I.
!
!    Output, integer I4_MODP, the nonnegative remainder when I is
!    divided by J.
!
  integer i
  integer i4_modp
  integer j

  if ( j == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'I4_MODP - Fatal error!'
    write ( *, '(a,i6)' ) '  I4_MODP ( I, J ) called with J = ', j
    stop
  end if

  i4_modp = mod ( i, j )

  if ( i4_modp < 0 ) then
    i4_modp = i4_modp + abs ( j )
  end if

  return
end
subroutine i4_swap ( i, j )

!*****************************************************************************80
!
!! I4_SWAP swaps two integer values.
!
!  Modified:
!
!    30 November 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer I, J.  On output, the values of I and
!    J have been interchanged.
!
  integer i
  integer j
  integer k

  k = i
  i = j
  j = k

  return
end
function i4_wrap ( ival, ilo, ihi )

!*****************************************************************************80
!
!! I4_WRAP forces an integer to lie between given limits by wrapping.
!
!  Example:
!
!    ILO = 4, IHI = 8
!
!    I  I4_WRAP
!
!    -2     8
!    -1     4
!     0     5
!     1     6
!     2     7
!     3     8
!     4     4
!     5     5
!     6     6
!     7     7
!     8     8
!     9     4
!    10     5
!    11     6
!    12     7
!    13     8
!    14     4
!
!  Modified:
!
!    15 July 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IVAL, an integer value.
!
!    Input, integer ILO, IHI, the desired bounds for the integer value.
!
!    Output, integer I4_WRAP, a "wrapped" version of IVAL.
!
  integer i4_modp
  integer i4_wrap
  integer ihi
  integer ilo
  integer ival
  integer wide

  wide = ihi + 1 - ilo

  if ( wide == 0 ) then
    i4_wrap = ilo
  else
    i4_wrap = ilo + i4_modp ( ival-ilo, wide )
  end if

  return
end
subroutine ihpsrt ( k, n, lda, a, map )

!*****************************************************************************80
!
!! IHPSRT uses heapsort on integer points in K-dimension.
!
!  Purpose:
!
!    Use heapsort to obtain the permutation of N K-dimensional
!    integer points so that the points are in lexicographic
!    increasing order.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer K, the dimension of points.
!
!    Input, integer N, the number of points.
!
!    Input, integer LDA, the leading dimension of array A in calling 
!    routine; K <= LDA.
!
!    Input, integer A(1:K,1:*), the array of at least N K-dimensional 
!    integer points.
!
!    Input/output, integer MAP(1:N), the points of A with indices MAP(1), 
!    MAP(2), ..., MAP(N) are to be sorted.
!    On output, elements are permuted so that A(*,MAP(1)) <=
!    A(*,MAP(2)) <= ... <= A(*,MAP(N)).
!
  integer lda
  integer n

  integer a(lda,*)
  integer i
  integer k
  integer map(n)
  integer t

  do i = n/2, 1, -1
    call isftdw ( i, n, k, lda, a, map )
  end do

  do i = n, 2, -1
    t = map(1)
    map(1) = map(i)
    map(i) = t
    call isftdw ( 1, i-1, k, lda, a, map )
  end do

  return
end
function iless ( k, p, q )

!*****************************************************************************80
!
!! ILESS determines whether a K-dimensional point P is lexically less than Q.
!
!  Purpose: 
!
!    Determine whether P is lexicographically less than Q.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer K, the dimension of points.
!
!    Input, integer P(1:K), Q(1:K), two points to compare.
!
!    Output, logical ILESS, is .TRUE. if P < Q, .FALSE. otherwise.
!
  integer k

  integer i
  logical iless
  integer p(k)
  integer q(k)

  do i = 1, k

    if ( p(i) /= q(i) ) then

      if ( p(i) < q(i) ) then
        iless = .true.
      else
        iless = .false.
      end if

      return

    end if

  end do

  iless = .false.

  return
end
subroutine i4mat_print ( lda, m, n, a, title )

!*****************************************************************************80
!
!! I4MAT_PRINT prints an I4MAT.
!
!  Modified:
!
!    08 July 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer LDA, the leading dimension of A.
!
!    Input, integer M, the number of rows in A.
!
!    Input, integer N, the number of columns in A.
!
!    Input, integer A(LDA,N), the matrix to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  integer lda
  integer n

  integer a(lda,n)
  integer i
  integer j
  integer jhi
  integer jlo
  integer m
  character ( len = * ) title

  if ( title /= ' ' ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  do jlo = 1, n, 10
    jhi = min ( jlo + 9, n )
    write ( *, '(a)' ) ' '
    write ( *, '(6x,10(i7))' ) ( j, j = jlo, jhi )
    write ( *, '(a)' ) ' '
    do i = 1, m
      write ( *, '(i6,10i7)' ) i, a(i,jlo:jhi)
    end do
  end do

  return
end
subroutine insed2 ( v, w, npolg, nvert, maxhv, maxpv, vcl, regnum, &
  hvl, pvl, iang, ierror )

!*****************************************************************************80
!
!! INSED2 inserts an edge into the head and polygon vertex lists.
!
!  Purpose: 
!
!    Insert the edge joining vertices V, W into the head vertex
!    list and polygon vertex list data structures.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer V, W, indices in PVL of vertices which are the endpoints
!    of an edge to be added to decomposition.
!
!    Input, integer NPOLG, the number of positions used in HVL array.
!
!    Input, integer NVERT, the number of positions used in PVL array.
!
!    Input, integer MAXHV, the maximum size available for HVL array.
!
!    Input, integer MAXPV, the maximum size available for PVL array.
!
!    Input, real ( kind = 8 ) VCL(1:2,1:*), the vertex coordinate list.
!
!    Input/output, integer REGNUM(1:NPOLG), the region numbers.
!
!    Input/output, integer HVL(1:NPOLG), the head vertex list.
!
!    Input/output, integer PVL(1:4,1:NVERT), real ( kind = 8 ) IANG(1:NVERT), 
!    the polygon vertex list and interior angles.
!
!    Output, integer IERROR, error flag.  On abnormal return,
!    IERROR is set to 4 or 5.
!
  integer maxhv
  integer maxpv

  real ( kind = 8 ) angle
  integer, parameter :: edgv = 4
  integer hvl(maxhv)
  integer i
  real ( kind = 8 ) iang(maxpv)
  integer ierror
  integer l
  integer, parameter :: loc = 1
  integer lv
  integer lw
  integer, parameter :: msglvl = 0
  integer npolg
  integer nvert
  integer, parameter :: polg = 2
  integer pvl(4,maxpv)
  integer regnum(maxhv)
  integer, parameter :: succ = 3
  integer v
  real ( kind = 8 ) vcl(2,*)
  integer vv
  integer w
  integer ww

  ierror = 0

  if ( maxhv <= npolg ) then
    ierror = 4
    return
  else if ( maxpv < nvert + 2 ) then
    ierror = 5
    return
  end if
!
!  Split linked list of vertices of the polygon containing vertices
!  V and W into two linked list of vertices of polygons with common
!  edge joining V and W.
!
  nvert = nvert + 2
  vv = nvert - 1
  ww = nvert
  lv = pvl(loc,v)
  lw = pvl(loc,w)
  pvl(loc,vv) = lv
  pvl(loc,ww) = lw
  pvl(polg,ww) = pvl(polg,v)
  pvl(succ,vv) = pvl(succ,v)
  pvl(succ,ww) = pvl(succ,w)
  pvl(succ,v) = ww
  pvl(succ,w) = vv
  pvl(edgv,vv) = pvl(edgv,v)
  pvl(edgv,ww) = pvl(edgv,w)
  pvl(edgv,v) = w
  pvl(edgv,w) = v

  if ( 0 < pvl(edgv,vv) ) then
    pvl(edgv,pvl(edgv,vv)) = vv
  end if

  if ( 0 < pvl(edgv,ww) ) then
    pvl(edgv,pvl(edgv,ww)) = ww
  end if

  l = pvl(loc,pvl(succ,vv))

  iang(vv) = angle ( vcl(1,lw), vcl(2,lw), vcl(1,lv), vcl(2,lv), vcl(1,l), &
    vcl(2,l) )

  iang(v) = iang(v) - iang(vv)
  l = pvl(loc,pvl(succ,ww))
  iang(ww) = angle ( vcl(1,lv), vcl(2,lv), vcl(1,lw), vcl(2,lw), vcl(1,l), &
    vcl(2,l) )

  iang(w) = iang(w) - iang(ww)
  npolg = npolg + 1
  i = vv

  do

    pvl(polg,i) = npolg
    i = pvl(succ,i)

    if ( i == vv ) then
      exit
    end if

  end do

  hvl(pvl(polg,v)) = v
  hvl(npolg) = vv
  regnum(npolg) = regnum(pvl(polg,v))

  if ( msglvl == 2 ) then
    write ( *, '(2i6,4f15.7)' ) v, w, vcl(1:2,lv), vcl(1:2,lw)
  end if

  return
end
subroutine insvr2 ( xi, yi, wp, nvc, nvert, maxvc, maxpv, vcl, pvl, &
  iang, w, ierror )

!*****************************************************************************80
!
!! INSVR2 inserts a point into the vertex coordinate and polygon vertex lists.
!
!  Purpose: 
!
!    Insert point (XI,YI) into vertex coordinate list and
!    polygon vertex list data structures.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) XI, YI, the coordinates of point to be inserted.
!
!    Input, integer WP, the index of vertex in PVL which is to be the
!    predecessor vertex of the inserted vertex.
!
!    Input/output, integer NVC, the number of positions used in VCL array.
!
!    Input/output, integer NVERT, the number of positions used in PVL array.
!
!    Input, integer MAXVC, the maximum size available for VCL array.
!
!    Input, integer MAXPV, the maximum size available for PVL array.
!
!    Input/output, real ( kind = 8 ) VCL(1:2,1:NVC), the vertex coordinate list.
!
!    Input/output, integer PVL(1:4,1:NVERT), real ( kind = 8 ) IANG(1:NVERT),
!    the polygon vertex list and interior angles.
!
!    Output, integer W, the index of inserted vertex in PVL.
!
!    Output, integer IERROR, error flag.  On abnormal return,
!    IERROR is set to 3 or 5.
!
  integer maxpv
  integer maxvc

  integer, parameter :: edgv = 4
  real ( kind = 8 ) iang(maxpv)
  integer ierror
  integer, parameter :: loc = 1
  integer nvc
  integer nvert
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  integer, parameter :: polg = 2
  integer pvl(4,maxpv)
  integer, parameter :: succ = 3
  real ( kind = 8 ) tol
  real ( kind = 8 ) vcl(2,maxvc)
  integer w
  integer wp
  integer ws
  integer ww
  integer wwp
  integer wws
  real ( kind = 8 ) xi
  real ( kind = 8 ) yi

  ierror = 0
  tol = 100.0D+00 * epsilon ( tol )

  if ( maxvc <= nvc ) then
    ierror = 3
    return
  end if

  if ( maxpv < nvert + 2 ) then
    ierror = 5
    return
  end if
!
!  Update linked list of vertices of polygon containing vertex WP.
!
  nvc = nvc + 1
  vcl(1,nvc) = xi
  vcl(2,nvc) = yi
  ws = pvl(succ,wp)
  nvert = nvert + 1
  w = nvert
  pvl(loc,w) = nvc
  pvl(polg,w) = pvl(polg,wp)
  pvl(succ,wp) = w
  pvl(succ,w) = ws
  iang(w) = pi
  pvl(edgv,w) = pvl(edgv,wp)
!
!  If edge containing (XI,YI) is shared by another polygon,
!  then also update linked list of vertices of that polygon.
!
  if ( 0 < pvl(edgv,wp) ) then
    wws = pvl(edgv,wp)
    wwp = pvl(succ,wws)
    nvert = nvert + 1
    ww = nvert
    pvl(loc,ww) = nvc
    pvl(polg,ww) = pvl(polg,wws)
    pvl(succ,wws) = ww
    pvl(succ,ww) = wwp
    iang(ww) = pi
    pvl(edgv,wp) = ww
    pvl(edgv,ww) = wp
    pvl(edgv,wws) = w
  end if

  return
end
subroutine intpg ( nvrt, xc, yc, ctrx, ctry, arpoly, hflag, umdf, wsq, nev, &
  ifv, listev, ivrt, edgval, vrtval, vcl, mdfint, mean, stdv, mdftr )

!*****************************************************************************80
!
!! INTPG integrates the mesh distribution function in a convex polygon.
!
!  Purpose: 
!
!    Compute integral of MDF2(X,Y) [heuristic mdf] or
!    UMDF(X,Y) [user-supplied mdf] in convex polygon.
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer NVRT, the number of vertices in polygon.
!
!    Input, real ( kind = 8 ) XC(0:NVRT), YC(0:NVRT), the coordinates of 
!    polygon vertices in counter clockwise order, translated so that 
!    centroid is at origin; (XC(0),YC(0)) = (XC(NVRT),YC(NVRT)).
!
!    Input, real ( kind = 8 ) CTRX, CTRY, the coordinates of centroid 
!    before translation.
!
!    Input, real ( kind = 8 ) ARPOLY, the area of polygon.
!
!    Input, logical HFLAG, is .TRUE. if heuristic mdf, .FALSE. if 
!    user-supplied mdf.
!
!    Input, external UMDF, the name of the user supplied MDF routine, of 
!    the form:
!
!      function umdf ( x, y )
!      double precision umdf
!      double precision x
!      double precision y
!
!    Input, real ( kind = 8 ) WSQ, the square of width of original polygon 
!    of decomposition.
!
!    Input, integer NEV, integer IFV, integer LISTEV(1:NEV), output from 
!    routine PRMDF2.
!
!    Input, IVRT(1:*), real ( kind = 8 ) EDGVAL(1:*), 
!    double precision VRTVAL(1:*), arrays output from DSMDF2;
!    if .NOT. HFLAG then only first array exists.
!
!    Input, real ( kind = 8 ) VCL(1:2,1:*), the vertex coordinate list.
!
!    Output, real ( kind = 8 ) MDFINT, the integral of MDF in polygon.
!
!    Output, real ( kind = 8 ) MEAN, the mean MDF value in polygon.
!
!    Output, real ( kind = 8 ) STDV, the standard deviation of MDF in polygon.
!
!    Output, real ( kind = 8 ) MDFTR(0:NVRT-1), the mean MDF value in each 
!    triangle of polygon;  triangles are determined by polygon vertices 
!    and centroid.
!
  integer nev
  integer, parameter :: nqpt = 3
  integer nvrt

  real ( kind = 8 ) areatr
  real ( kind = 8 ) arpoly
  real ( kind = 8 ) ctrx
  real ( kind = 8 ) ctry
  real ( kind = 8 ) d
  real ( kind = 8 ) edgval(*)
  logical hflag
  integer i
  integer ifv
  integer ivrt(*)
  integer j
  integer k
  integer kp1
  integer l
  integer listev(nev)
  integer m
  real ( kind = 8 ) mdfint
  real ( kind = 8 ) mdfsqi
  real ( kind = 8 ) mdftr(0:nvrt-1)
  real ( kind = 8 ) mean
  real ( kind = 8 ), save, dimension ( 3, nqpt ) :: qc = reshape ( (/ &
    0.6666666666666666D+00, 0.1666666666666667D+00, 0.1666666666666667D+00, &
    0.1666666666666667D+00, 0.6666666666666666D+00, 0.1666666666666667D+00, &
    0.1666666666666667D+00, 0.1666666666666667D+00, 0.6666666666666666D+00/), &
    (/ 3, nqpt /) )
  real ( kind = 8 ) s
  real ( kind = 8 ) stdv
  real ( kind = 8 ) sum1
  real ( kind = 8 ) sum2
  real ( kind = 8 ) temp
  real ( kind = 8 ) umdf
  real ( kind = 8 ) val
  real ( kind = 8 ) vcl(2,*)
  real ( kind = 8 ) vrtval(*)
  real ( kind = 8 ) wsq
  real ( kind = 8 ), save, dimension ( nqpt ) :: wt = &
    (/ 0.3333333333333333D+00, 0.3333333333333333D+00, 0.3333333333333333D+00 /)
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) x1
  real ( kind = 8 ) xc(0:nvrt)
  real ( kind = 8 ) xx
  real ( kind = 8 ) y
  real ( kind = 8 ) y0
  real ( kind = 8 ) y1
  real ( kind = 8 ) yc(0:nvrt)
  real ( kind = 8 ) yy

  external umdf
!
!  NQPT is number of quadrature points for numerical integration in triangle.
!  WT(I) is weight of Ith quadrature point.
!  QC(1:3,I) are barycentric coordinates of Ith quadrature point.
!
  mdfint = 0.0D+00
  mdfsqi = 0.0D+00

  do l = 0, nvrt-1

    areatr = 0.5D+00 * ( xc(l) * yc(l+1) - xc(l+1) * yc(l) )
    sum1 = 0.0D+00
    sum2 = 0.0D+00

    do m = 1, nqpt

      xx = qc(1,m) * xc(l) + qc(2,m) * xc(l+1)
      yy = qc(1,m) * yc(l) + qc(2,m) * yc(l+1)
      if ( hflag ) then
!
!             VAL = MDF2(XX+CTRX,YY+CTRY,WSQ,NEV,IFV,LISTEV,IVRT, &
!               EDGVAL,VRTVAL,VCL)
!
!  Insert code for function MDF2 to reduce number of calls.
!
         x = xx + ctrx
         y = yy + ctry
           s = wsq
           do i = 1, nev
            k = listev(i)
            if ( k < 0 ) then
               k = -k
               d = ( vcl(1,k) - x )**2 + ( vcl(2,k) - y )**2
               d = max ( 0.25D+00 * d, vrtval(k) )
               s = min ( s, d )
            else
               kp1 = k + 1
               if ( i == nev .and. 0 < ifv ) then
                 kp1 = ifv
               end if
               j = ivrt(kp1)
               x0 = x - vcl(1,j)
               y0 = y - vcl(2,j)
               x1 = vcl(1,ivrt(k)) - vcl(1,j)
               y1 = vcl(2,ivrt(k)) - vcl(2,j)

               if ( x0 * x1 + y0 * y1 <= 0.0D+00 ) then
                 d = x0**2 + y0**2
               else
                 x0 = x0 - x1
                 y0 = y0 - y1
                 if ( 0.0D+00 <= x0 * x1 + y0 * y1 ) then
                   d = x0**2 + y0**2
                 else
                   d = ( x1 * y0 - y1 * x0 )**2 / ( x1**2 + y1**2 )
                 end if
               end if

               d = max ( 0.25D+00 * d, edgval(k) )
               s = min ( s, d )
            end if
           end do
           val = 1.0D+00 / s
      else
        val = umdf ( xx+ctrx, yy+ctry )
      end if

      temp = wt(m) * val
      sum1 = sum1 + temp
      sum2 = sum2 + temp * val

    end do

    mdftr(l) = sum1
    mdfint = mdfint + sum1 * areatr
    mdfsqi = mdfsqi + sum2 * areatr

  end do

  mean = mdfint / arpoly
  stdv = mdfsqi / arpoly - mean**2
  stdv = sqrt ( max ( stdv, 0.0D+00 ) )

  return
end
subroutine inttri ( nvrt, xc, yc, h, ibot, costh, sinth, ldv, nvc, ntri,  &
  maxvc, maxti, maxcw, vcl, til, ncw, cwalk, ierror )

!*****************************************************************************80
!
!! INTTRI generates triangles inside a convex polygon.
!
!  Purpose: 
!
!    Generate triangles inside convex polygon using quasi-uniform grid of
!    spacing H.  It is assumed that the diameter of the polygon is parallel 
!    to the Y axis.
!
!  Modified:
!
!    02 May 2001
!
!  Author:
!
!    Barry Joe, 
!    Department of Computing Science, 
!    University of Alberta,
!    Edmonton, Alberta, Canada  T6G 2H1
!
!  Reference:
!
!    Barry Joe,
!    GEOMPACK - a software package for the generation of meshes
!    using geometric algorithms,
!    Advances in Engineering Software,
!    Volume 13, pages 325-331, 1991.
!
!  Parameters:
!
!    Input, integer NVRT, the number of vertices on the boundary of
!    convex polygon.
!
!    Input, real ( kind = 8 ) XC(0:NVRT), YC(0:NVRT), the vertex coordinates 
!    in counter clockwise order; (XC(0),YC(0)) = (XC(NVRT),YC(NVRT)).
!
!    Input, real ( kind = 8 ) H, the spacing of mesh vertices in polygon.
!
!    Input, integer IBOT, the index of bottom vertex; diameter contains vertices
!    (XC(0),YC(0)) and (XC(IBOT),YC(IBOT)).
!
!    Input, real ( kind = 8 ) COSTH, SINTH; COS(THETA), SIN(THETA) where 
!    THETA in [-PI,PI] is rotation angle to get diameter parallel to y-axis.
!
!    Input, integer LDV, the leading dimension of VCL in calling routine.
!
!    Input/output, integer NVC, the number of coordinates or positions 
!    used in VCL array.
!
!    Input/output, integer NTRI, the number of triangles or positions 
!    used in TIL.
!
!    Input, integer MAXVC, the maximum size available for VCL array.
!
!    Input, integer MAXTI, the maximum size available for TIL array.
!
!    Input, integer MAXCW, the maximum size available for CWALK array; 
!    assumed to be at least 6*(1 + INT((YC(0) - YC(IBOT))/H)).
!
!    Input/output, real ( kind = 8 ) VCL(1:2,1:NVC), the vertex coordinate list.
!
!    Input/output, integer TIL(1:3,1:NTRI), the triangle incidence list.
!
!    Output, integer NCW, the number of mesh vertices in closed walk, 
!    except NCW = 0 for 1 vertex.
!
!    Output, integer CWALK(0:NCW), indices in VCL of mesh vertices of closed
!    walk; CWALK(0) = CWALK(NCW)
!
!    Output, integer IERROR, error flag.  On abnormal return,
!    IERROR is set to 3, 9, or 10
!
  integer ldv
  integer maxcw
  integer maxti
  integer maxvc
  integer nvrt

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) costh
  integer cwalk(0:maxcw)
  real ( kind = 8 ) cy
  real ( kind = 8 ) h
  integer i
  integer ibot
  integer ierror
  integer il
  integer im1l
  integer im1r
  integer ir
  integer j
  integer k
  integer l
  integer l0
  integer l1
  integer lw
  integer m
  integer n
  integer ncw
  integer ntri
  integer nvc
  integer p
  integer r
  integer r0
  integer r1
  integer rw
  real ( kind = 8 ) sinth
  real ( kind = 8 ) sy
  integer til(3,maxti)
  real ( kind = 8 ) tol
  real ( kind = 8 ) vcl(ldv,maxvc)
  real ( kind = 8 ) x
  real ( kind = 8 ) xc(0:nvrt)
  real ( kind = 8 ) xj
  real ( kind = 8 ) xk
  real ( kind = 8 ) xl
  real ( kind = 8 ) xm1l
  real ( kind = 8 ) xm1r
  real ( kind = 8 ) xr
  real ( kind = 8 ) y
  real ( kind = 8 ) yc(0:nvrt)

  ierror = 0
  tol = 100.0D+00 * epsilon ( tol )

  n = int ( ( yc(0) - yc(ibot) ) / h )
  y = yc(0) - 0.5D+00 * ( yc(0) - yc(ibot ) - real ( n, kind = 8 ) * h )
  l = 0
  r = nvrt

  do i = 0, n
!
!  Determine left and right x-coordinates of polygon for
!  scan line with y-coordinate Y, and generate mesh vertices.
!
    do while ( y < yc(l+1) )
      l = l + 1
    end do

    do while ( y < yc(r-1) )
      r = r - 1
    end do

    xl = xc(l) + ( xc(l+1) - xc(l) ) * ( y - yc(l) ) / ( yc(l+1) - yc(l) )
    xr = xc(r) + ( xc(r-1) - xc(r) ) * ( y - yc(r) ) / ( yc(r-1) - yc(r) )
    m = int ( ( xr - xl ) / h )
    x = xl + 0.5D+00 * ( xr - xl - real ( m, kind = 8 ) * h )

    if ( maxvc < nvc + m + 1 ) then
      ierror = 3
      return
    end if

    cy = costh * y
    sy = sinth * y
    il = nvc + 1
    xl = x

    do j = 0, m
      nvc = nvc + 1
      vcl(1,nvc) = costh * x + sy
      vcl(2,nvc) = cy - sinth * x
      x = x + h
    end do

    ir = nvc
    xr = x - h

    if ( n == 0 ) then

      ncw = 0
      cwalk(0) = nvc
      return

    else if ( i == 0 ) then

      lw = 0
      cwalk(lw) = il
      rw = maxcw + 1

      do j = il, ir
        rw = rw - 1
        cwalk(rw) = j
      end do

      go to 100

    end if
!
!  Generate triangles between scan lines Y+H and Y.
!
    a = max ( xl, xm1l )
    b = min ( xr, xm1r )

    if ( xm1l == a ) then
      l0 = im1l
      x = ( xm1l - xl ) / h
      j = int(x + tol)
      if ( abs ( x - real ( j, kind = 8 ) ) <= tol ) then
        j = j - 1
      end if
      if ( j < 0 ) then
        j = 0
      end if
      l1 = il + j
    else
      l1 = il
      x = ( xl - xm1l ) / h
      j = int ( x + tol )
      if ( abs ( x - real ( j, kind = 8 ) ) <= tol ) then
        j = j - 1
      end if
      if ( j < 0 ) then
        j = 0
      end if
      l0 = im1l + j
    end if

    if ( xm1r == b ) then
      r0 = im1r
      x = ( xr - xm1r ) / h
      j = int ( x + tol )
      if ( abs ( x - real ( j, kind = 8 ) ) <= tol ) then
        j = j - 1
      end if
      if ( j < 0 ) then
        j = 0
      end if
      r1 = ir - j
    else
      r1 = ir
      x = ( xm1r - xr ) / h
      j = int ( x + tol )
      if ( abs ( x - real ( j, kind = 8 ) ) <= tol ) then
        j = j - 1
      end if
      if ( j < 0 ) then
        j = 0
      end if
      r0 = im1r - j
    end if

    if ( l0 < r0 .or. l1 < r1 ) then

      j = l0
      k = l1
      xj = xm1l + real ( j-im1l, kind = 8 ) * h
      xk = xl + real ( k - il, kind = 8 ) * h

      do

        if ( k < r1 .and. ( xk <= xj .or. j == r0 ) ) then
          p = k
          k = k + 1
          xk = xk + h
        else
          p = j
          j = j + 1
          xj = xj + h
        end if

        ntri = ntri + 1

        if ( maxti < ntri ) then
          ierror = 9
          return
        end if

        til(1,ntri) = j
        til(2,ntri) = p
        til(3,ntri) = k

        if ( r0 <= j .and. r1 <= k ) then
          exit
        end if

      end do

    end if
!
!  Generate paths of closed walk between scan lines Y+H and Y.
!
    if ( xm1l < xl ) then
      do j = im1l+1, l0
        lw = lw + 1
        cwalk(lw) =