function angle_deg_2d ( x1, y1, x2, y2, x3, y3 )
!
!*******************************************************************************
!
!! ANGLE_DEG_2D returns the angle swept out between two rays in 2D.
!
!
!  Discussion:
!
!    Except for the zero angle case, it should be true that
!
!      ANGLE_DEG_2D(X1,Y1,X2,Y2,X3,Y3)
!    + ANGLE_DEG_2D(X3,Y3,X2,Y2,X1,Y1) = 360.0E+00
!
!  Modified:
!
!    14 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, X3, Y3, define the rays
!    ( X1-X2, Y1-Y2 ) and ( X3-X2, Y3-Y2 ) which in turn define the
!    angle, counterclockwise from ( X1-X2, Y1-Y2 ).
!
!    Output, real ANGLE_DEG_2D, the angle swept out by the rays, measured
!    in degrees.  0 <= ANGLE_DEG_2D < 360.  If either ray has zero length,
!    then ANGLE_DEG_2D is set to 0.
!
  implicit none
!
  real angle_deg_2d
  real angle_rad_2d
  real pi
  real radians_to_degrees
  real x
  real x1
  real x2
  real x3
  real y
  real y1
  real y2
  real y3
!
  x = ( x1 - x2 ) * ( x3 - x2 ) + ( y1 - y2 ) * ( y3 - y2 )
  y = ( x1 - x2 ) * ( y3 - y2 ) - ( y1 - y2 ) * ( x3 - x2 )

  if ( x == 0.0E+00 .and. y == 0.0E+00 ) then

    angle_deg_2d = 0.0E+00

  else

    angle_rad_2d = atan2 ( y, x )

    if ( angle_rad_2d < 0.0E+00 ) then
      angle_rad_2d = angle_rad_2d + 2.0E+00 * pi ( )
    end if

    angle_deg_2d = radians_to_degrees ( angle_rad_2d )

  end if

  return
end
function angle_rad_2d ( x1, y1, x2, y2, x3, y3 )
!
!*******************************************************************************
!
!! ANGLE_RAD_2D returns the angle swept out between two rays in 2D.
!
!
!  Discussion:
!
!    Except for the zero angle case, it should be true that
!
!      ANGLE_RAD_2D(X1,Y1,X2,Y2,X3,Y3)
!    + ANGLE_RAD_2D(X3,Y3,X2,Y2,X1,Y1) = 2 * PI
!
!  Modified:
!
!    19 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, X3, Y3, define the rays
!    ( X1-X2, Y1-Y2 ) and ( X3-X2, Y3-Y2 ) which in turn define the
!    angle, counterclockwise from ( X1-X2, Y1-Y2 ).
!
!    Output, real ANGLE_RAD_2D, the angle swept out by the rays,
!    in radians.  0 <= ANGLE_RAD_2D < 2 PI.  If either ray has zero
!    length, then ANGLE_RAD_2D is set to 0.
!
  implicit none
!
  real angle_rad_2d
  real pi
  real x
  real x1
  real x2
  real x3
  real y
  real y1
  real y2
  real y3
!
  x = ( x1 - x2 ) * ( x3 - x2 ) + ( y1 - y2 ) * ( y3 - y2 )
  y = ( x1 - x2 ) * ( y3 - y2 ) - ( y1 - y2 ) * ( x3 - x2 )

  if ( x == 0.0E+00 .and. y == 0.0E+00 ) then

    angle_rad_2d = 0.0E+00

  else

    angle_rad_2d = atan2 ( y, x )

    if ( angle_rad_2d < 0.0E+00 ) then
      angle_rad_2d = angle_rad_2d + 2.0E+00 * pi ( )
    end if

  end if

  return
end
subroutine circle_dia2imp_2d ( x1, y1, x2, y2, r, cx, cy )
!
!*******************************************************************************
!
!! CIRCLE_DIA2IMP_2D converts a diameter to an implicit circle in 2D.
!
!
!  Discussion:
!
!    The diameter form of a circle is:
!
!      (X1,Y1) and (X2,Y2) are endpoints of a diameter.
!
!    The implicit form of a circle in 2D is:
!
!      (X-CX)**2 + (Y-CY)**2 = R**2
!
!  Modified:
!
!    27 January 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, are the X and Y coordinates
!    of two points which form a diameter of the circle.
!
!    Output, real R, the computed radius of the circle.
!
!    Output, real CX, CY, the computed center of the circle.
!
  implicit none
!
  real enorm0_2d
  real r
  real x1
  real x2
  real cx
  real y1
  real y2
  real cy
!
  r = 0.5E+00 * enorm0_2d ( x1, y1, x2, y2 )

  cx = 0.5E+00 * ( x1 + x2 )
  cy = 0.5E+00 * ( y1 + y2 )

  return
end
subroutine circle_exp2imp_2d ( x1, y1, x2, y2, x3, y3, r, cx, cy )
!
!*******************************************************************************
!
!! CIRCLE_EXP2IMP_2D converts a circle from explicit to implicit form in 2D.
!
!
!  Formula:
!
!    The explicit form of a circle in 2D is:
!
!      The circle passing through (X1,Y1), (X2,Y2), (X3,Y3).
!
!    The implicit form of a circle in 2D is:
!
!      (X-CX)**2 + (Y-CY)**2 = R**2
!
!  Discussion:
!
!    Any three points define a circle, as long as they don't lie on a straight
!    line.  (If the points do lie on a straight line, we could stretch the
!    definition of a circle to allow an infinite radius and a center at
!    some infinite point.)
!
!    Instead of the formulas used here, you can use the linear system
!    approach in the routine TRIANGLE_OUTCIRCLE_2D.
!
!    The diameter of the circle can be found by solving a 2 by 2 linear system.
!    This is because the vectors P2 - P1 and P3 - P1 are secants of the circle,
!    and each forms a right triangle with the diameter.  Hence, the dot product
!    of P2 - P1 with the diameter is equal to the square of the length
!    of P2 - P1, and similarly for P3 - P1.  These two equations determine the
!    diameter vector originating at P1.
!
!  Reference:
!
!    Joseph O'Rourke,
!    Computational Geometry,
!    Cambridge University Press,
!    Second Edition, 1998, page 187.
!
!  Modified:
!
!    12 June 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, X3, Y3, are the X and Y coordinates
!    of three points that lie on the circle.  These points should be
!    distinct, and not collinear.
!
!    Output, real R, the radius of the circle.  Normally, R will be positive.
!    R is returned as -1 in the unlikely event that the points are
!    numerically collinear.
!
!    Output, real CX, CY, the center of the circle.
!
  implicit none
!
  real a
  real b
  real c
  real d
  real e
  real f
  real g
  real r
  real x1
  real x2
  real x3
  real cx
  real y1
  real y2
  real y3
  real cy
!
  a = x2 - x1
  b = y2 - y1
  c = x3 - x1
  d = y3 - y1

  e = a * ( x1 + x2 ) + b * ( y1 + y2 )
  f = c * ( x1 + x3 ) + d * ( y1 + y3 )
!
!  Our formula is:
!
!    G = a * ( d - b ) - b * ( c - a )
!
!  but we get slightly better results using the original data.
!
  g = a * ( y3 - y2 ) - b * ( x3 - x2 )
!
!  We check for collinearity.  A more useful check would compare the
!  absolute value of G to a small quantity.
!
  if ( g == 0.0E+00 ) then
    cx = 0.0E+00
    cy = 0.0E+00
    r = -1.0E+00
    return
  end if
!
!  The center is halfway along the diameter vector from (X1,Y1).
!
  cx = 0.5E+00 * ( d * e - b * f ) / g
  cy = 0.5E+00 * ( a * f - c * e ) / g
!
!  Knowing the center, the radius is now easy to compute.
!
  r = sqrt ( ( x1 - cx )**2 + ( y1 - cy )**2 )

  return
end
subroutine circle_imp_contains_point_2d ( r, cx, cy, x, y, inside )
!
!*******************************************************************************
!
!! CIRCLE_IMP_CONTAINS_POINT_2D determines if an implicit circle contains a point in 2D.
!
!
!  Formula:
!
!    An implicit circle in 2D satisfies the equation:
!
!      ( X - CX )**2 + ( Y - CY )**2 = R**2
!
!  Modified:
!
!    21 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real R, the radius of the circle.
!
!    Input, real CX, CY, the coordinates of the center of the circle.
!
!    Input, real X, Y, the point to be checked.
!
!    Output, logical INSIDE, is TRUE if the point is inside or on the circle,
!    FALSE otherwise.
!
  implicit none
!
  real enormsq0_2d
  logical inside
  real r
  real x
  real cx
  real y
  real cy
!
  if ( enormsq0_2d ( x, y, cx, cy ) <= r * r ) then
    inside = .true.
  else
    inside = .false.
  end if

  return
end
function cross0_2d ( x0, y0, x1, y1, x2, y2 )
!
!*******************************************************************************
!
!! CROSS0_2D finds the cross product of (P1-P0) and (P2-P0) in 2D.
!
!
!  Discussion:
!
!    Strictly speaking, the vectors lie in the (X,Y) plane, and
!    the cross product here is a vector in the Z direction.
!
!    The vectors are specified with respect to a basis point P0.
!    We are computing the normal to the triangle (P0,P1,P2).
!
!  Modified:
!
!    19 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real X3, Y3, Z3, the cross product (P1-P0) x (P2-P0), or
!    (X1-X0,Y1-Y0,Z1-Z0) x (X2-X0,Y2-Y0,Z2-Z0).
!
!    Input, real X0, Y0, X1, Y1, X2, Y2, the coordinates of the three
!    points.  The basis point P0 is (X0,Y0).
!
!    Output, real CROSS0_2D, the Z component of the cross product
!    (X1-X0,Y1-Y0,0) x (X2-X0,Y2-Y0,0).
!
  implicit none
!
  real cross0_2d
  real x0
  real x1
  real x2
  real y0
  real y1
  real y2
!
  cross0_2d = ( x1 - x0 ) * ( y2 - y0 ) - ( y1 - y0 ) * ( x2 - x0 )

  return
end
function enorm0_2d ( x0, y0, x1, y1 )
!
!*******************************************************************************
!
!! ENORM0_2D computes the Euclidean norm of (P1-P0) in 2D.
!
!
!  Modified:
!
!    16 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X0, Y0, X1, Y1, the coordinates of the points P0 and P1.
!
!    Output, real ENORM0_2D, the Euclidean norm of (P1-P0).
!
  implicit none
!
  real enorm0_2d
  real x0
  real x1
  real y0
  real y1
!
  enorm0_2d = sqrt ( ( x1 - x0 )**2 + ( y1 - y0 )**2 )

  return
end
function enormsq0_2d ( x0, y0, x1, y1 )
!
!*******************************************************************************
!
!! ENORMSQ0_2D computes the square of the Euclidean norm of (P1-P0) in 2D.
!
!
!  Modified:
!
!    06 February 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X0, Y0, X1, Y1, the coordinates of the points
!    P0 and P1.
!
!    Output, real ENORMSQ0_2D, the square of the Euclidean norm of (P1-P0).
!
  implicit none
!
  real enormsq0_2d
  real x0
  real x1
  real y0
  real y1
!
  enormsq0_2d = ( x1 - x0 )**2 + ( y1 - y0 )**2

  return
end
function i_modp ( i, j )
!
!*******************************************************************************
!
!! I_MODP returns the nonnegative remainder of integer division.
!
!
!  Formula:
!
!    If
!      NREM = I_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, I_MODP(A,360) is between 0 and 360, always.
!
!  Examples:
!
!        I     J     MOD  I_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 I_MODP, the nonnegative remainder when I is
!    divided by J.
!
  implicit none
!
  integer i
  integer i_modp
  integer j
!
  if ( j == 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'I_MODP - Fatal error!'
    write ( *, * ) '  I_MODP ( I, J ) called with J = ', j
    stop
  end if

  i_modp = mod ( i, j )

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

  return
end
subroutine i_random ( ilo, ihi, i )
!
!*******************************************************************************
!
!! I_RANDOM returns a random integer in a given range.
!
!
!  Modified:
!
!    23 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ILO, IHI, the minimum and maximum acceptable values.
!
!    Output, integer I, the randomly chosen integer.
!
  implicit none
!
  logical, save :: seed = .false.
  integer i
  integer ihi
  integer ilo
  real r
  real rhi
  real rlo
!
  if ( .not. seed ) then
    call random_seed
    seed = .true.
  end if
!
!  Pick a random number in (0,1).
!
  call random_number ( harvest = r )
!
!  Set a real interval [RLO,RHI] which contains the integers [ILO,IHI],
!  each with a "neighborhood" of width 1.
!
  rlo = real ( ilo ) - 0.5E+00
  rhi = real ( ihi ) + 0.5E+00
!
!  Set I to the integer that is nearest the scaled value of R.
!
  i = nint ( ( 1.0E+00 - r ) * rlo + r * rhi )
!
!  In case of oddball events at the boundary, enforce the limits.
!
  i = max ( i, ilo )
  i = min ( i, ihi )

  return
end
subroutine i_swap ( i, j )
!
!*******************************************************************************
!
!! I_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.
!
  implicit none
!
  integer i
  integer j
  integer k
!
  k = i
  i = j
  j = k

  return
end
function i_wrap ( ival, ilo, ihi )
!
!*******************************************************************************
!
!! I_WRAP forces an integer to lie between given limits by wrapping.
!
!
!  Example:
!
!    ILO = 4, IHI = 8
!
!    I  I_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 I_WRAP, a "wrapped" version of IVAL.
!
  implicit none
!
  integer i_modp
  integer i_wrap
  integer ihi
  integer ilo
  integer ival
  integer wide
!
  wide = ihi + 1 - ilo

  if ( wide == 0 ) then
    i_wrap = ilo
  else
    i_wrap = ilo + i_modp ( ival-ilo, wide )
  end if

  return
end
subroutine ij_next ( i, j, n )
!
!*******************************************************************************
!
!! IJ_NEXT returns the next matrix index.
!
!
!  Discussion:
!
!    For N = 3, the sequence of indices returned is:
!
!      (1,1), (1,2), (1,3), (2,1), (2,2), (2,3), (3,1), (3,2), (3,3), (0,0).
!
!    Note that once the value (N,N) is returned, the next value returned
!    will be (0,0).
!
!  Modified:
!
!    28 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer I, J.  On input, the current pair of indices.
!    On output, the next pair of indices.  If either index is illegal on
!    input, the output value of (I,J) will be (1,1).
!
!    Input, integer N, the maximum value for I and J.
!
  implicit none
!
  integer i
  integer j
  integer n
!
  if ( n < 1 ) then
    i = 0
    j = 0
    return
  end if

  if ( i < 1 .or. i > n .or. j < 1 .or. j > n ) then
    i = 1
    j = 1
    return
  end if

  if ( j < n ) then
    j = j + 1
  else if ( i < n ) then
    i = i + 1
    j = 1
  else
    i = 0
    j = 0
  end if

  return
end
subroutine ij_next_gt ( i, j, n )
!
!*******************************************************************************
!
!! IJ_NEXT_GT returns the next matrix index, with the constraint that I < J.
!
!
!  Discussion:
!
!    For N = 5, the sequence of indices returned is:
!
!      (1,2), (1,3), (1,4), (1,5), (2,3), (2,4), (2,5), (3,4), (3,5), (4,5).
!
!  Modified:
!
!    28 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer I, J.  On input, the current pair of indices.
!    On output, the next pair of indices.  If either index is illegal on
!    input, the output value of (I,J) will be (1,2).
!
!    Input, integer N, the maximum value for I and J.
!    A value of N less than 2 is nonsense.
!
  implicit none
!
  integer i
  integer j
  integer n
!
  if ( n < 2 ) then
    i = 0
    j = 0
    return
  end if

  if ( i < 1 .or. i > n .or. j < 1 .or. j > n .or. j <= i ) then
    i = 1
    j = 2
    return
  end if

  if ( j < n ) then
    j = j + 1
  else if ( i < n - 1 ) then
    i = i + 1
    j = i + 1
  else
    i = 0
    j = 0
  end if

  return
end
subroutine imat_print ( lda, m, n, a, title )
!
!*******************************************************************************
!
!! IMAT_PRINT prints an integer matrix.
!
!
!  Modified:
!
!    08 May 2000
!
!  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.
!
  implicit none
!
  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 ( *, * ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

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

  return
end
subroutine ivec_frac ( n, a, k, iafrac )
!
!*******************************************************************************
!
!! IVEC_FRAC searches for the K-th smallest element in an N-vector.
!
!
!  Discussion:
!
!    Hoare's algorithm is used.
!
!  Modified:
!
!    17 July 2000
!
!  Parameters:
!
!    Input, integer N, the number of elements of A.
!
!    Input/output, integer A(N), array to search.  On output,
!    the elements of A have been somewhat rearranged.
!
!    Input, integer K, the fractile to be sought.  If K = 1, the
!    minimum entry is sought.  If K = N, the maximum is sought.
!    Other values of K search for the entry which is K-th in size.
!    K must be at least 1, and no greater than N.
!
!    Output, integer IAFRAC, the value of the K-th fractile of A.
!
  implicit none
!
  integer n
!
  integer i
  integer a(n)
  integer iafrac
  integer iryt
  integer iw
  integer ix
  integer j
  integer k
  integer left
!
  if ( n <= 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'IVEC_FRAC  - Fatal error!'
    write ( *, * ) '  Illegal nonpositive value of N = ', n
    stop
  end if

  if ( k <= 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'IVEC_FRAC  - Fatal error!'
    write ( *, * ) '  Illegal nonpositive value of K = ', k
    stop
  end if

  if ( k > n ) then
    write ( *, * ) ' '
    write ( *, * ) 'IVEC_FRAC  - Fatal error!'
    write ( *, * ) '  Illegal K > N, K = ', k
    stop
  end if

  left = 1
  iryt = n

  do

    if ( left >= iryt ) then
      iafrac = a(k)
      exit
    end if

    ix = a(k)
    i = left
    j = iryt

    do

      if ( i > j ) then

        if ( j < k ) then
          left = i
        end if

        if ( k < i ) then
          iryt = j
        end if

        exit

      end if
!
!  Find I so that IX <= A(I).
!
      do while ( a(i) < ix )
        i = i + 1
      end do
!
!  Find J so that A(J) <= IX.
!
      do while ( ix < a(j) )
        j = j - 1
      end do

      if ( i <= j ) then
        call i_swap ( a(i), a(j) )
        i = i + 1
        j = j - 1
      end if

    end do

  end do

  return
end
subroutine ivec_heap_a ( n, a )
!
!*******************************************************************************
!
!! IVEC_HEAP_A reorders an array of integers into an ascending heap.
!
!
!  Definition:
!
!    An ascending heap is an array A with the property that, for every index J,
!    A(J) <= A(2*J) and A(J) <= A(2*J+1), (as long as the indices
!    2*J and 2*J+1 are legal).
!
!  Diagram:
!
!                  A(1)
!                /      \
!            A(2)         A(3)
!          /     \        /  \
!      A(4)       A(5)  A(6) A(7)
!      /  \       /   \
!    A(8) A(9) A(10) A(11)
!
!  Reference:
!
!    A Nijenhuis and H Wilf,
!    Combinatorial Algorithms,
!    Academic Press, 1978, second edition,
!    ISBN 0-12-519260-6.
!
!  Modified:
!
!    20 March 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the size of the input array.
!
!    Input/output, integer A(N).
!    On input, an unsorted array.
!    On output, the array has been reordered into a heap.
!
  implicit none
!
  integer n
!
  integer a(n)
  integer i
  integer ifree
  integer key
  integer m
!
!  Only nodes N/2 down to 1 can be "parent" nodes.
!
  do i = n/2, 1, -1
!
!  Copy the value out of the parent node.
!  Position IFREE is now "open".
!
    key = a(i)
    ifree = i

    do
!
!  Positions 2*IFREE and 2*IFREE + 1 are the descendants of position
!  IFREE.  (One or both may not exist because they exceed N.)
!
      m = 2 * ifree
!
!  Does the first position exist?
!
      if ( m > n ) then
        exit
      end if
!
!  Does the second position exist?
!
      if ( m + 1 <= n ) then
!
!  If both positions exist, take the smaller of the two values,
!  and update M if necessary.
!
        if ( a(m+1) < a(m) ) then
          m = m + 1
        end if

      end if
!
!  If the small descendant is smaller than KEY, move it up,
!  and update IFREE, the location of the free position, and
!  consider the descendants of THIS position.
!
      if ( a(m) >= key ) then
        exit
      end if

      a(ifree) = a(m)
      ifree = m

    end do
!
!  Once there is no more shifting to do, KEY moves into the free spot.
!
    a(ifree) = key

  end do

  return
end
subroutine ivec_heap_d ( n, a )
!
!*******************************************************************************
!
!! IVEC_HEAP_D reorders an array of integers into an descending heap.
!
!
!  Definition:
!
!    A descending heap is an array A with the property that, for every index J,
!    A(J) >= A(2*J) and A(J) >= A(2*J+1), (as long as the indices
!    2*J and 2*J+1 are legal).
!
!  Diagram:
!
!                  A(1)
!                /      \
!            A(2)         A(3)
!          /     \        /  \
!      A(4)       A(5)  A(6) A(7)
!      /  \       /   \
!    A(8) A(9) A(10) A(11)
!
!  Reference:
!
!    A Nijenhuis and H Wilf,
!    Combinatorial Algorithms,
!    Academic Press, 1978, second edition,
!    ISBN 0-12-519260-6.
!
!  Modified:
!
!    15 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the size of the input array.
!
!    Input/output, integer A(N).
!    On input, an unsorted array.
!    On output, the array has been reordered into a heap.
!
  implicit none
!
  integer n
!
  integer a(n)
  integer i
  integer ifree
  integer key
  integer m
!
!  Only nodes N/2 down to 1 can be "parent" nodes.
!
  do i = n/2, 1, -1
!
!  Copy the value out of the parent node.
!  Position IFREE is now "open".
!
    key = a(i)
    ifree = i

    do
!
!  Positions 2*IFREE and 2*IFREE + 1 are the descendants of position
!  IFREE.  (One or both may not exist because they exceed N.)
!
      m = 2 * ifree
!
!  Does the first position exist?
!
      if ( m > n ) then
        exit
      end if
!
!  Does the second position exist?
!
      if ( m + 1 <= n ) then
!
!  If both positions exist, take the larger of the two values,
!  and update M if necessary.
!
        if ( a(m+1) > a(m) ) then
          m = m + 1
        end if

      end if
!
!  If the large descendant is larger than KEY, move it up,
!  and update IFREE, the location of the free position, and
!  consider the descendants of THIS position.
!
      if ( a(m) <= key ) then
        exit
      end if

      a(ifree) = a(m)
      ifree = m

    end do
!
!  Once there is no more shifting to do, KEY moves into the free spot IFREE.
!
    a(ifree) = key

  end do

  return
end
subroutine ivec_heap_d_extract ( n, a, val )
!
!*******************************************************************************
!
!! IVEC_HEAP_D_EXTRACT extracts the maximum value from a descending heap.
!
!
!  Discussion:
!
!    In other words, the routine finds the maximum value in the
!    heap, returns that value to the user, deletes that value from
!    the heap, and restores the heap to its proper form.
!
!    This is one of three functions needed to model a priority queue.
!
!  Reference:
!
!    Thomas Cormen, Charles Leiserson, Ronald Rivest,
!    Introduction to Algorithms,
!    MIT Press, page 150.
!
!  Modified:
!
!    28 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer N, the number of items in the heap.
!
!    Input/output, integer A(N), the heap.
!
!    Output, integer VAL, the item of maximum value, which has been
!    removed from the heap.
!
  implicit none
!
  integer a(*)
  integer n
  integer val
!
  if ( n < 1 ) then
    write ( *, * ) ' '
    write ( *, * ) 'I_HEAP_D_EXTRACT - Fatal error!'
    write ( *, * ) '  The heap is empty.'
    stop
  end if
!
!  Get the maximum value.
!
  val = a(1)

  if ( n == 1 ) then
    n = 0
    return
  end if
!
!  Shift the last value down.
!
  a(1) = a(n)
!
!  Restore the heap structure.
!
  n = n - 1
  call ivec_sort_heap_d ( n, a )

  return
end
subroutine ivec_heap_d_insert ( n, a, val )
!
!*******************************************************************************
!
!! IVEC_HEAP_D_INSERT inserts a new value into a descending heap.
!
!
!  Discussion:
!
!    This is one of three functions needed to model a priority queue.
!
!  Reference:
!
!    Thomas Cormen, Charles Leiserson, Ronald Rivest,
!    Introduction to Algorithms,
!    MIT Press, page 150.
!
!  Modified:
!
!    28 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer N, the number of items in the heap.
!
!    Input/output, integer A(N), the heap.
!
!    Input, integer VAL, the value to be inserted.
!
  implicit none
!
  integer a(*)
  integer i
  integer n
  integer parent
  integer val
!
  n = n + 1
  i = n

  do while ( i > 1 )

    parent = i / 2

    if ( a(parent) >= val ) then
      exit
    end if

    a(i) = a(parent)
    i = parent

  end do

  a(i) = val

  return
end
subroutine ivec_heap_d_max ( n, a, val_max )
!
!*******************************************************************************
!
!! IVEC_HEAP_D_MAX returns the maximum value in a descending heap of integers.
!
!
!  Discussion:
!
!    This is one of three functions needed to model a priority queue.
!
!  Reference:
!
!    Thomas Cormen, Charles Leiserson, Ronald Rivest,
!    Introduction to Algorithms,
!    MIT Press, page 150.
!
!  Modified:
!
!    28 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of items in the heap.
!
!    Input, integer A(N), the heap.
!
!    Output, integer VAL_MAX, the maximum value in the heap.
!
  implicit none
!
  integer n
!
  integer a(n)
  integer val_max
!
  val_max = a(1)

  return
end
subroutine ivec_identity ( n, a )
!
!*******************************************************************************
!
!! IVEC_IDENTITY sets an integer vector to the identity vector A(I)=I.
!
!
!  Modified:
!
!    09 November 2000
!
!  Parameters:
!
!    Input, integer N, the number of elements of A.
!
!    Output, integer A(N), the array to be initialized.
!
  implicit none
!
  integer n
!
  integer a(n)
  integer i
!
  do i = 1, n
    a(i) = i
  end do

  return
end
subroutine ivec_median ( n, a, median )
!
!*******************************************************************************
!
!! IVEC_MEDIAN returns the median of an unsorted integer vector.
!
!
!  Discussion:
!
!    Hoare's algorithm is used.  The values of the vector are
!    rearranged by this routine.
!
!  Modified:
!
!    18 September 2000
!
!  Parameters:
!
!    Input, integer N, the number of elements of A.
!
!    Input/output, integer A(N), the array to search.  On output,
!    the order of the elements of A has been somewhat changed.
!
!    Output, integer MEDIAN, the value of the median of A.
!
  implicit none
!
  integer n
!
  integer a(n)
  integer k
  integer median
!
  k = ( n + 1 ) / 2

  call ivec_frac ( n, a, k, median )

  return
end
subroutine ivec_pop ( n, x, stack1_max, stack1_num, stack1, stack2_max, &
  stack2_num, stack2 )
!
!*******************************************************************************
!
!! IVEC_POP pops an integer vector off of a stack.
!
!
!  Discussion:
!
!    If there are no more objects in the stack, N is returned as -1.
!
!  Modified:
!
!    09 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer N, the dimension of the vector.
!
!    Output, integer X(*), the value of the vector.
!
!    Input, integer STACK1_MAX, the maximum size of STACK1.
!
!    Input/output, integer STACK1_NUM, the current size of STACK1.
!
!    Input/output, integer STACK1(STACK1_MAX), the vector dimension stack.
!
!    Input, integer STACK2_MAX, the maximum size of STACK2.
!
!    Input/output, integer STACK2_NUM, the current size of STACK2.
!
!    Input/output, integer STACK2(STACK2_MAX), the vector value stack.
!
  implicit none
!
  integer n
  integer stack1_max
  integer stack2_max
!
  integer stack1(stack1_max)
  integer stack1_num
  integer stack2(stack2_max)
  integer stack2_num
  integer x(*)
!
  if ( stack1_num < 1 ) then
    n = -1
    return
  end if

  n = stack1(stack1_num)
  stack1_num = stack1_num - 1

  stack2_num = stack2_num - n
  x(1:n) = stack2(stack2_num+1:stack2_num+n)

  return
end
subroutine ivec_print ( n, a, title )
!
!*******************************************************************************
!
!! IVEC_PRINT prints an integer vector.
!
!
!  Modified:
!
!    16 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, integer A(N), the vector to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none
!
  integer n
!
  integer a(n)
  integer i
  character ( len = * ) title
!
  if ( title /= ' ' ) then
    write ( *, * ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, * ) ' '
  do i = 1, n
    write ( *, '(i6,i10)' ) i, a(i)
  end do

  return
end
subroutine ivec_push ( n, x, stack1_max, stack1_num, stack1, stack2_max, &
  stack2_num, stack2 )
!
!*******************************************************************************
!
!! IVEC_PUSH pushes an integer vector onto a stack.
!
!
!  Discussion:
!
!    STACK1 contains a list of the dimensions of the objects stored.
!    Therefore, STACK1_MAX should be at least as big as the maximum number
!    of objects to be considered.
!
!    STACK2 contains the values of the objects.  Therefore, STACK2_MAX
!    should probably be as big as the maximum total length of the maximum
!    number of objects stored.
!
!    On first call, the user should have set STACK1_NUM and STACK2_NUM to zero.
!
!  Modified:
!
!    09 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the dimension of the vector.  N may be zero.
!
!    Input, integer X(N), the value of the vector.
!
!    Input, integer STACK1_MAX, the maximum size of STACK1.
!
!    Input/output, integer STACK1_NUM, the current size of STACK1.
!
!    Input/output, integer STACK1(STACK1_MAX), the vector dimension stack.
!
!    Input, integer STACK2_MAX, the maximum size of STACK2.
!
!    Input/output, integer STACK2_NUM, the current size of STACK2.
!
!    Input/output, integer STACK2(STACK2_MAX), the vector value stack.
!
  implicit none
!
  integer n
  integer stack1_max
  integer stack2_max
!
  integer stack1(stack1_max)
  integer stack1_num
  integer stack2(stack2_max)
  integer stack2_num
  integer x(n)
!
  if ( n < 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'IVEC_PUSH - Fatal error!'
    write ( *, * ) '  Input dimension N is negative.'
    stop
  end if

  if ( stack1_num + 1 > stack1_max ) then
    write ( *, * ) ' '
    write ( *, * ) 'IVEC_PUSH - Fatal error!'
    write ( *, * ) '  Exceeding size of stack #1.'
    stop
  end if

  if ( stack2_num + n > stack2_max ) then
    write ( *, * ) ' '
    write ( *, * ) 'IVEC_PUSH - Fatal error!'
    write ( *, * ) '  Exceeding size of stack #2.'
    stop
  end if

  stack1_num = stack1_num + 1
  stack1(stack1_num) = n

  stack2(stack2_num+1:stack2_num+n) = x(1:n)
  stack2_num = stack2_num + n

  return
end
subroutine ivec_reverse ( n, a )
!
!*******************************************************************************
!
!! IVEC_REVERSE reverses the elements of an integer vector.
!
!
!  Example:
!
!    Input:
!
!      N = 5,
!      A = ( 11, 12, 13, 14, 15 ).
!
!    Output:
!
!      A = ( 15, 14, 13, 12, 11 ).
!
!  Modified:
!
!    26 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the array.
!
!    Input/output, integer A(N), the array to be reversed.
!
  implicit none
!
  integer n
!
  integer a(n)
  integer i
!
  do i = 1, n/2
    call i_swap ( a(i), a(n+1-i) )
  end do

  return
end
subroutine ivec_sort_heap_d ( n, a )
!
!*******************************************************************************
!
!! IVEC_SORT_HEAP_D descending sorts an integer array using heap sort.
!
!
!  Reference:
!
!    A Nijenhuis and H Wilf,
!    Combinatorial Algorithms,
!    Academic Press, 1978, second edition,
!    ISBN 0-12-519260-6.
!
!  Modified:
!
!    15 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the array.
!
!    Input/output, integer A(N).
!    On input, the array to be sorted;
!    On output, the array has been sorted.
!
  implicit none
!
  integer n
!
  integer a(n)
  integer n1
!
  if ( n <= 1 ) then
    return
  end if
!
!  1: Put A into ascending heap form.
!
  call ivec_heap_a ( n, a )
!
!  2: Sort A.
!
!  The smallest object in the heap is in A(1).
!  Move it to position A(N).
!
  call i_swap ( a(1), a(n) )
!
!  Consider the diminished heap of size N1.
!
  do n1 = n-1, 2, -1
!
!  Restore the heap structure of A(1) through A(N1).
!
    call ivec_heap_a ( n1, a )
!
!  Take the smallest object from A(1) and move it to A(N1).
!
    call i_swap ( a(1), a(n1) )

  end do

  return
end
subroutine ivec_split_unsort ( n, a, split, isplit )
!
!*******************************************************************************
!
!! IVEC_SPLIT_UNSORT "splits" an unsorted vector based on a splitting value.
!
!
!  Discussion:
!
!    If the vector is already sorted, it is simpler to do a binary search
!    on the data than to call this routine.
!
!    The vector is not assumed to be sorted before input, and is not
!    sorted during processing.  If sorting is not needed, then it is
!    more efficient to use this routine.
!
!  Modified:
!
!    18 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of elements of A.
!
!    Input/output, integer A(N), the array to split.  On output,
!    all the entries of A that are less than or equal to SPLIT
!    are in A(1:ISPLIT).
!
!    Input, integer SPLIT, the value used to split the vector.
!    It is not necessary that any value of A actually equal SPLIT.
!
!    Output, integer ISPLIT, indicates the position of the last
!    entry of the split vector that is less than or equal to SPLIT.
!
  implicit none
!
  integer n
!
  integer a(n)
  integer i
  integer i1
  integer i2
  integer i3
  integer isplit
  integer j1
  integer j2
  integer j3
  integer split
!
!  Partition the vector into A1, A2, A3, where
!    A1 = A(I1:J1) holds values <= SPLIT,
!    A2 = A(I2:J2) holds untested values,
!    A3 = A(I3:J3) holds values > SPLIT.
!
  i1 = 1
  j1 = 0

  i2 = 1
  j2 = n

  i3 = n+1
  j3 = n
!
!  Pick the next item from A2, and move it into A1 or A3.
!  Adjust indices appropriately.
!
  do i = 1, n

    if ( a(i2) <= split ) then
      i2 = i2 + 1
      j1 = j1 + 1
    else
      call i_swap ( a(i2), a(i3-1) )
      i3 = i3 - 1
      j2 = j2 - 1
    end if

  end do

  isplit = j1

  return
end
subroutine line_exp2imp_2d ( x1, y1, x2, y2, a, b, c )
!
!*******************************************************************************
!
!! LINE_EXP2IMP_2D converts an explicit line to implicit form in 2D.
!
!
!  Formula:
!
!    The explicit form of a line in 2D is:
!
!      (X1,Y1), (X2,Y2).
!
!    The implicit form of a line in 2D is:
!
!      A * X + B * Y + C = 0
!
!  Modified:
!
!    24 January 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2.  (X1,Y1) and (X2,Y2) are
!    two points on the line. (X1,Y1) must be different
!    from (X2,Y2).
!
!    Output, real A, B, C, three coefficients which describe
!    the line that passes through (X1,Y1) and (X2,Y2).
!
  implicit none
!
  real a
  real b
  real c
  real x1
  real x2
  real y1
  real y2
!
!  Take care of degenerate cases.
!
  if ( x1 == x2 .and. y1 == y2 ) then
    write ( *, * ) ' '
    write ( *, * ) 'LINE_EXP2IMP_2D - Fatal error!'
    write ( *, * ) '  (X1,Y1) = (X2,Y2)'
    write ( *, * ) '  (X1,Y1) = ', x1, y1
    write ( *, * ) '  (X2,Y2) = ', x2, y2
    stop
  end if

  a = y2 - y1
  b = x1 - x2
  c = x2 * y1 - x1 * y2

  return
end
subroutine line_exp_point_dist_2d ( x1, y1, x2, y2, x, y, dist )
!
!*******************************************************************************
!
!! LINE_EXP_POINT_DIST_2D: distance ( explicit line, point ) in 2D.
!
!
!  Formula:
!
!    The explicit form of a line in 2D is:
!
!      (X1,Y1), (X2,Y2).
!
!  Modified:
!
!    27 January 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2.  (X1,Y1) and (X2,Y2) are two points on
!    the line.
!
!    Input, real X, Y, the point whose distance from the line is
!    to be measured.
!
!    Output, real DIST, the distance from the point to the line.
!
  implicit none
!
  real bot
  real dist
  real dot
  real enorm0_2d
  real enormsq0_2d
  real t
  real x
  real xn
  real x1
  real x2
  real y
  real yn
  real y1
  real y2
!
  bot = enormsq0_2d ( x1, y1, x2, y2 )

  if ( bot == 0.0E+00 ) then

    xn = x1
    yn = y1
!
!  (P-P1) dot (P2-P1) = Norm(P-P1) * Norm(P2-P1) * Cos(Theta).
!
!  (P-P1) dot (P2-P1) / Norm(P-P1)**2 = normalized coordinate T
!  of the projection of (P-P1) onto (P2-P1).
!
  else

    dot = ( x - x1 ) * ( x2 - x1 ) + ( y - y1 ) * ( y2 - y1 )

    t = dot / bot

    xn = x1 + t * ( x2 - x1 )
    yn = y1 + t * ( y2 - y1 )

  end if

  dist = enorm0_2d ( xn, yn, x, y )

  return
end
subroutine line_exp_point_dist_signed_2d ( x1, y1, x2, y2, x, y, dist_signed )
!
!*******************************************************************************
!
!! LINE_EXP_POINT_DIST_SIGNED_2D: signed distance ( explicit line, point ) in 2D.
!
!
!  Discussion:
!
!    The signed distance has two interesting properties:
!
!    *  The absolute value of the signed distance is the
!        usual (Euclidean) distance.
!
!    *  Points with signed distance 0 lie on the line,
!       points with a negative signed distance lie on one side
!         of the line,
!       points with a positive signed distance lie on the
!         other side of the line.
!
!    Assuming that C is nonnegative, then if a point is a positive
!    distance away from the line, it is on the same side of the
!    line as the point (0,0), and if it is a negative distance
!    from the line, it is on the opposite side from (0,0).
!
!  Modified:
!
!    06 February 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, define the two points
!    (X1,Y1) and (X2,Y2) that determine the line.
!
!    Input, real X, Y, the point (X,Y) whose signed distance is desired.
!
!    Output, real DIST_SIGNED, the signed distance from the point to the line.
!
  implicit none
!
  real a
  real b
  real c
  real dist_signed
  real x
  real x1
  real x2
  real y
  real y1
  real y2
!
!  Convert the line to implicit form.
!
  call line_exp2imp_2d ( x1, y1, x2, y2, a, b, c )
!
!  Compute the signed distance from the point to the line.
!
  dist_signed = ( a * x + b * y + c ) / sqrt ( a * a + b * b )

  return
end
subroutine line_seg_contains_point_2d ( x1, y1, x2, y2, x3, y3, u, v )
!
!*******************************************************************************
!
!! LINE_SEG_CONTAINS_POINT_2D reports if a line segment contains a point in 2D.
!
!
!  Discussion:
!
!    In exact arithmetic, point P3 = (X3,Y3) is on the line segment between
!    P1=(X1,Y1) and P2=(X2,Y2) if and only if 0 <= U <= 1 and V = 0.
!
!  Modified:
!
!    28 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, the endpoints P1 and P2 of a line segment.
!
!    Input, real X3, Y3, a point P3 to be tested.
!
!    Output, real U, the coordinate of (X3,Y3) along the axis from
!    with origin at P1 and unit at P2.
!
!    Output, real V, the magnitude of the off-axis portion of the
!    vector P3-P1, measured in units of (P2-P1).
!
  implicit none
!
  real u
  real unit
  real v
  real x1
  real x2
  real x3
  real y1
  real y2
  real y3
!
  unit = sqrt ( ( x2 - x1 )**2 + ( y2 - y1 )**2 )

  if ( unit == 0.0E+00 ) then

    if ( x3 == x1 .and. y3 == y1 ) then
      u = 0.5E+00
      v = 0.0E+00
    else
      u = 0.5E+00
      v = huge ( 1.0E+00 )
    end if

  else

    u = ( ( x3 - x1 ) * ( x2 - x1 ) + ( y3 - y1 ) * ( y2 - y1 ) ) / unit**2

    v = sqrt ( ( ( u - 1.0E+00 ) * x1 - u * x2 + x3 )**2 &
             + ( ( u - 1.0E+00 ) * y1 - u * y2 + y3 )**2 ) / unit

  end if

  return
end
subroutine line_seg_vec_int_2d ( n, x1, y1, x2, y2, i, j, flag, xint, yint )
!
!*******************************************************************************
!
!! LINE_SEG_VEC_INT_2D computes intersections of a set of line segments.
!
!
!  Discussion:
!
!    This is an implementation of the relatively inefficient algorithm
!    for computing all intersections of a set of line segments.
!
!    This is an "incremental" code, which returns as soon as it finds
!    a single intersection.  To find the next intersection, simply call
!    again.
!
!  Modified:
!
!    28 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of line segments.
!
!    Input, real X1(N), Y1(N), X2(N), Y2(N), the coordinates of the
!    endpoints of the line segments.
!
!    Input/output, integer I, J, used to keep track of the computation.
!    On first call with a given set of line segments, set I = J = 0.
!    On return with FLAG = 1, I and J record the indices of the
!    line segments whose intersection has just been found.  To find the
!    next intersection, simply call again, but do not alter I and J.
!
!    Output, integer FLAG:
!    0, no more intersections, the computation is done.
!    1, an intersection was detected between segments I and J.
!
!    Output, real XINT, YINT, the location of an intersection of line
!    segments I and J, if FLAG is 1.
!
  implicit none
!
  integer n
!
  integer flag
  integer i
  integer j
  real x1(n)
  real x2(n)
  real xint
  real y1(n)
  real y2(n)
  real yint
!
  do

    call ij_next_gt ( i, j, n )

    if ( i == 0 ) then
      flag = 0
      exit
    end if

    call lines_seg_int_2d ( x1(i), y1(i), x2(i), y2(i), x1(j), y1(j), &
      x2(j), y2(j), flag, xint, yint )

    if ( flag == 1 ) then
      return
    end if

  end do

  return
end
subroutine lines_exp_int_2d ( ival, x, y, x1, y1, x2, y2, x3, y3, x4, y4 )
!
!*******************************************************************************
!
!! LINES_EXP_INT_2D determines where two explicit lines intersect in 2D.
!
!
!  Formula:
!
!    The explicit form of a line in 2D is:
!
!      (X1,Y1), (X2,Y2).
!
!  Modified:
!
!    28 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IVAL, reports on the intersection.
!
!     0, no intersection, the lines may be parallel or degenerate.
!     1, one intersection point, returned in X, Y.
!     2, infinitely many intersections, the lines are identical.
!
!    Output, real X, Y, if IVAl = 1, then X, Y contains
!    the intersection point.  Otherwise, X = 0, Y = 0.
!
!    Input, real X1, Y1, X2, Y2, define the first line.
!
!    Input, real X3, Y3, X4, Y4, define the second line.
!
  implicit none
!
  real a1
  real a2
  real b1
  real b2
  real c1
  real c2
  integer ival
  logical point_1
  logical point_2
  real x
  real x1
  real x2
  real x3
  real x4
  real y
  real y1
  real y2
  real y3
  real y4
!
  ival = 0
  x = 0.0E+00
  y = 0.0E+00
!
!  Check whether either line is a point.
!
  if ( x1 == x2 .and. y1 == y2 ) then
    point_1 = .true.
  else
    point_1 = .false.
  end if

  if ( x3 == x4 .and. y3 == y4 ) then
    point_2 = .true.
  else
    point_2 = .false.
  end if
!
!  Convert the lines to ABC format.
!
  if ( .not. point_1 ) then
    call line_exp2imp_2d ( x1, y1, x2, y2, a1, b1, c1 )
  end if

  if ( .not. point_2 ) then
    call line_exp2imp_2d ( x3, y3, x4, y4, a2, b2, c2 )
  end if
!
!  Search for intersection of the lines.
!
  if ( point_1 .and. point_2 ) then
    if ( x1 == x3 .and. y1 == y3 ) then
      ival = 1
      x = x1
      y = y1
    end if
  else if ( point_1 ) then
    if ( a2 * x1 + b2 * y1 == c2 ) then
      ival = 1
      x = x1
      y = y1
    end if
  else if ( point_2 ) then
    if ( a1 * x3 + b1 * y3 == c1 ) then
      ival = 1
      x = x3
      y = y3
    end if
  else
    call lines_imp_int_2d ( a1, b1, c1, a2, b2, c2, ival, x, y )
  end if

  return
end
subroutine lines_imp_int_2d ( a1, b1, c1, a2, b2, c2, ival, x, y )
!
!*******************************************************************************
!
!! LINES_IMP_INT_2D determines where two implicit lines intersect in 2D.
!
!
!  Formula:
!
!    The implicit form of a line in 2D is:
!
!      A * X + B * Y + C = 0
!
!  Modified:
!
!    16 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real A1, B1, C1, define the first line.
!    At least one of A1 and B1 must be nonzero.
!
!    Input, real A2, B2, C2, define the second line.
!    At least one of A2 and B2 must be nonzero.
!
!    Output, integer IVAL, reports on the intersection.
!
!    -1, both A1 and B1 were zero.
!    -2, both A2 and B2 were zero.
!     0, no intersection, the lines are parallel.
!     1, one intersection point, returned in X, Y.
!     2, infinitely many intersections, the lines are identical.
!
!    Output, real X, Y, if IVAL = 1, then X, Y contains
!    the intersection point.  Otherwise, X = 0, Y = 0.
!
  implicit none
!
  real a(2,2)
  real a1
  real a2
  real b(2,2)
  real b1
  real b2
  real c1
  real c2
  real det
  integer ival
  real x
  real y
!
  x = 0.0E+00
  y = 0.0E+00
!
!  Refuse to handle degenerate lines.
!
  if ( a1 == 0.0E+00 .and. b1 == 0.0E+00 ) then
    ival = -1
    return
  else if ( a2 == 0.0E+00 .and. b2 == 0.0E+00 ) then
    ival = -2
    return
  end if
!
!  Set up a linear system, and compute its inverse.
!
  a(1,1) = a1
  a(1,2) = b1
  a(2,1) = a2
  a(2,2) = b2

  call rmat2_inverse ( a, b, det )
!
!  If the inverse exists, then the lines intersect.
!  Multiply the inverse times -C to get the intersection point.
!
  if ( det /= 0.0E+00 ) then

    ival = 1
    x = - b(1,1) * c1 - b(1,2) * c2
    y = - b(2,1) * c1 - b(2,2) * c2
!
!  If the inverse does not exist, then the lines are parallel
!  or coincident.  Check for parallelism by seeing if the
!  C entries are in the same ratio as the A or B entries.
!
  else

    ival = 0

    if ( a1 == 0.0E+00 ) then
      if ( b2 * c1 == c2 * b1 ) then
        ival = 2
      end if
    else
      if ( a2 * c1 == c2 * a1 ) then
        ival = 2
      end if
    end if

  end if

  return
end
subroutine lines_seg_dist_2d ( x1, y1, x2, y2, x3, y3, x4, y4, flag, x5, y5 )
!
!*******************************************************************************
!
!! LINES_SEG_DIST_2D computes the distance of two line segments in 2D.
!
!
!  Modified:
!
!    28 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, the endpoints of the first segment.
!
!    Input, real X3, Y3, X4, Y4, the endpoints of the second segment.
!
!    Output, integer FLAG, records the results.
!    0, the line segments do not intersect.
!    1, the line segments intersect.
!
!    Output, real X5, Y5.
!    If FLAG = 0, X5 = Y5 = 0.
!    If FLAG = 1, then (X5,Y5) is a point of intersection.
!
  implicit none
!
  integer flag
  integer ival
  real u
  real v
  real x1
  real x2
  real x3
  real x4
  real x5
  real y1
  real y2
  real y3
  real y4
  real y5
!
  x5 = 0.0E+00
  y5 = 0.0E+00
!
!  Find the intersection of the two lines.
!
  call lines_exp_int_2d ( ival, x5, y5, x1, y1, x2, y2, x3, y3, x4, y4 )

  if ( ival == 0 ) then
    flag = 0
    return
  end if
!
!  Is the point on the first segment?
!
  call line_seg_contains_point_2d ( x1, y1, x2, y2, x5, y5, u, v )

  if ( u < 0 .or. 1.0E+00 < u .or. v > 0.001E+00 ) then
    flag = 0
    return
  end if
!
!  Is the point on the second segment?
!
  call line_seg_contains_point_2d ( x3, y3, x4, y4, x5, y5, u, v )

  if ( u < 0 .or. 1.0E+00 < u .or. v > 0.001E+00 ) then
    flag = 0
    return
  end if

  flag = 1

  return
end
subroutine lines_seg_int_1d ( x1, x2, x3, x4, flag, x5, x6 )
!
!*******************************************************************************
!
!! LINES_SEG_INT_1D computes the intersection of two line segments in 1D.
!
!
!  Modified:
!
!    07 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, X2, the endpoints of the first segment.
!
!    Input, real X3, X4, the endpoints of the second segment.
!
!    Output, integer FLAG, records the results.
!    0, the line segments do not intersect.
!    1, the line segments intersect.
!
!    Output, real X5, X6, the endpoints of the intersection segment.
!    If FLAG = 0, X5 = X6 = 0.
!
  implicit none
!
  integer flag
  real x1
  real x2
  real x3
  real x4
  real x5
  real x6
  real y1
  real y2
  real y3
  real y4
!
  y1 = min ( x1, x2 )
  y2 = max ( x1, x2 )
  y3 = min ( x3, x4 )
  y4 = max ( x3, x4 )
!
  flag = 0
  x5 = 0.0E+00
  x6 = 0.0E+00

  if ( y4 < y1 ) then
    return
  else if ( y3 > y2 ) then
    return
  end if

  flag = 1
  x5 = max ( y1, y3 )
  x6 = min ( y2, y4 )

  return
end
subroutine lines_seg_int_2d ( x1, y1, x2, y2, x3, y3, x4, y4, flag, x5, y5 )
!
!*******************************************************************************
!
!! LINES_SEG_INT_2D computes the intersection of two line segments in 2D.
!
!
!  Modified:
!
!    28 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, the endpoints of the first segment.
!
!    Input, real X3, Y3, X4, Y4, the endpoints of the second segment.
!
!    Output, integer FLAG, records the results.
!    0, the line segments do not intersect.
!    1, the line segments intersect.
!
!    Output, real X5, Y5.
!    If FLAG = 0, X5 = Y5 = 0.
!    If FLAG = 1, then (X5,Y5) is a point of intersection.
!
  implicit none
!
  integer flag
  integer ival
  real u
  real v
  real x1
  real x2
  real x3
  real x4
  real x5
  real y1
  real y2
  real y3
  real y4
  real y5
!
  x5 = 0.0E+00
  y5 = 0.0E+00
!
!  Find the intersection of the two lines.
!
  call lines_exp_int_2d ( ival, x5, y5, x1, y1, x2, y2, x3, y3, x4, y4 )

  if ( ival == 0 ) then
    flag = 0
    return
  end if
!
!  Is the point on the first segment?
!
  call line_seg_contains_point_2d ( x1, y1, x2, y2, x5, y5, u, v )

  if ( u < 0 .or. 1.0E+00 < u .or. v > 0.001E+00 ) then
    flag = 0
    return
  end if
!
!  Is the point on the second segment?
!
  call line_seg_contains_point_2d ( x3, y3, x4, y4, x5, y5, u, v )

  if ( u < 0 .or. 1.0E+00 < u .or. v > 0.001E+00 ) then
    flag = 0
    return
  end if

  flag = 1

  return
end
subroutine perm_print ( n, p, title )
!
!*******************************************************************************
!
!! PERM_PRINT prints a permutation.
!
!
!  Example:
!
!    Input:
!
!      P = 7 2 4 1 5 3 6
!
!    Printed output:
!
!      "This is the permutation:"
!
!      1 2 3 4 5 6 7
!      7 2 4 1 5 3 6
!
!  Modified:
!
!    25 July 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of objects permuted.
!
!    Input, integer P(N), the permutation, in standard index form.
!
!    Input, character ( len = * ) TITLE, an optional title.
!    If no title is supplied, then only the permutation is printed.
!
  implicit none
!
  integer, parameter :: inc = 20
!
  integer n
!
  integer i
  integer ihi
  integer ilo
  integer p(n)
  character ( len = * ) title
!
  if ( len_trim ( title ) /= 0 ) then

    write ( *, * ) ' '
    write ( *, '(a)' ) trim ( title )

    do ilo = 1, n, inc
      ihi = min ( n, ilo + inc - 1 )
      write ( *, * ) ' '
      write ( *, '(20i4)' ) ( i, i = ilo, ihi )
      write ( *, '(20i4)' ) p(ilo:ihi)
    end do

  else

    do ilo = 1, n, inc
      ihi = min ( n, ilo + inc - 1 )
      write ( *, '(20i4)' ) p(ilo:ihi)
    end do

  end if

  return
end
subroutine perm_random ( n, p )
!
!*******************************************************************************
!
!! PERM_RANDOM returns a random permutation.
!
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, page 78.
!
!  Modified:
!
!    12 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of items to permute.
!
!    Output, integer P(N), a permutation of the numbers from 1 to N.
!
  implicit none
!
  integer n
!
  integer i
  integer j
  integer p(n)
!
  call ivec_identity ( n, p )

  do i = n, 2, -1
    call i_random ( 1, i, j )
    call i_swap ( p(i), p(j) )
  end do

  return
end
function pi ( )
!
!*******************************************************************************
!
!! PI returns the value of pi.
!
!
!  Modified:
!
!    04 December 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real PI, the value of pi.
!
  implicit none
!
  real pi
!
  pi = 3.14159265358979323846264338327950288419716939937510E+00

  return
end
subroutine points_convex_hull_cubic_2d ( nu, ux, uy, nv, v )
!
!*******************************************************************************
!
!! POINTS_CONVEX_HULL_CUBIC_2D computes the convex hull of 2D points.
!
!
!  Discussion:
!
!    The algorithm used requires time that is cubic in the number of input
!    points.
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, pages 2-8.
!
!  Modified:
!
!    30 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NU, the number of points.
!
!    Input, real UX(NU), UY(NU), the coordinates of the points.
!
!    Output, integer NV, the number of vertices in the convex hull.
!
!    Output, integer V(NU), the NV vertices that form the convex hull,
!    in counter clockwise order.   (The arrays should be dimensioned to NU.  
!    The value of NV is guaranteed to be no greater than NU.)
!
  implicit none
!
  integer nu
!
  integer b(nu)
  real cross
  real cross0_2d
  integer e(nu)
  integer i
  integer j
  integer k
  integer match
  integer nv
  real ux(nu)
  real uy(nu)
  logical valid
  integer v(nu)
!
  nv = 0

  do i = 1, nu

    do j = 1, nu
      if ( i /= j ) then
        valid = .true.
        do k = 1, nu
          if ( k /= i .and. k /= j ) then
!
!  Compute the cross product: P(J)-P(I) x P(K)-P(I)
!  to determine if P(K) is strictly to the left of the line from P(I) to P(J).
!
!  I'll have to think if CROSS should be positive or negative...
!
            cross = cross0_2d ( ux(i), uy(i), ux(j), uy(j), ux(k), uy(k) )

            if ( cross > 0.0E+00 ) then
              valid = .false.
              exit
            end if

          end if

        end do
!
!  Add the line from P(I) to P(J) to the list.
!
        if ( valid ) then
          nv = nv + 1
          b(nv) = i
          e(nv) = j
        end if

      end if

    end do

  end do
!
!  From the unordered edge list of vertex pairs,
!  construct an ordered list of vertices.
!
  v(1) = b(1)
  match = e(1)

  do k = 2, nv

    v(k) = 0

    do j = k, nv
      if ( b(j) == match ) then
        v(k) = b(j)
        match = e(j)
        call i_swap ( b(j), b(k) )
        call i_swap ( e(j), e(k) )
        exit
      end if
    end do

    if ( v(k) == 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'POINTS_2D_CONVEX_HULL_CUBIC - Fatal error!'
      write ( *, * ) '  Algorithm failed  for K = ', k
      write ( *, * ) ' '
      write ( *, * ) '  B, E arrays:'
      write ( *, * ) ' '
      do i = 1, nv
        write ( *, '(2i6)' ) b(i), e(i)
      end do
      stop
    end if

  end do

  if ( match /= v(1) ) then
    write ( *, * ) ' '
    write ( *, * ) 'POINTS_2D_CONVEX_HULL_CUBIC - Fatal error!'
    write ( *, * ) '  Algorithm failed to link last node to first.'
    stop
  end if

!
!  Reverse the order.
!
  call ivec_reverse ( nv, v )
 
  return
end
subroutine points_convex_hull_nlogh_2d ( nu, ux, uy, nv, v )
!
!*******************************************************************************
!
!! POINTS_CONVEX_HULL_NLOGH_2D computes the convex hull of 2D points.
!
!
!  Discussion:
!
!    The work involved is N*log(H), where N is the number of points, and H is
!    the number of points that are on the hull.
!
!  Modified:
!
!    10 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NU, the number of nodes.
!
!    Input, real UX(NU), UY(NU), the X and Y coordinates of the nodes.
!
!    Output, integer NV, the number of nodes that lie on the convex hull.
!
!    Output, integer V(NU).  Entries 1 through NV contain
!    the indices of the nodes that form the convex hull, in order.
!
  implicit none
!
  integer nu
!
  real ai
  real angle_deg_2d
  real angmax
  real di
  real dr
  real enorm0_2d
  integer i
  integer iq
  integer ir
  integer istart
  integer nv
  real ux(nu)
  real uy(nu)
  integer v(nu)
  real xp
  real xq
  real xr
  real yp
  real yq
  real yr
!
  if ( nu < 1 ) then
    nv = 0
    return
  end if
!
!  If NU = 1, the hull is the point.
!
  if ( nu == 1 ) then
    nv = 1
    v(1) = 1
    return
  end if
!
!  If NU = 2, then the convex hull is either the two distinct points,
!  or possibly a single (repeated) point.
!
  if ( nu == 2 ) then

    if ( ux(1) /= ux(2) .or. uy(1) /= uy(2) ) then
      nv = 2
      v(1) = 1
      v(2) = 2
    else
      nv = 1
      v(1) = 1
    end if

    return

  end if
!
!  Find the leftmost point, and take the bottom-most in a tie.
!  Call it "Q".
!
  iq = 1
  do i = 2, nu
    if ( ux(i) < ux(iq) .or. &
       ( ux(i) == ux(iq) .and. uy(i) < uy(iq) ) ) then
      iq = i
    end if
  end do

  xq = ux(iq)
  yq = uy(iq)
!
!  Remember the starting point.
!
  istart = iq
  nv = 1
  v(1) = iq
!
!  For the first point, make a dummy previous point, 1 unit south,
!  and call it "P".
!
  xp = xq
  yp = yq - 1.0E+00
!
!  Now, having old point P, and current point Q, find the new point R
!  so the angle PQR is maximal.
!
!  Watch out for the possibility that the two nodes are identical.
!
  do

    ir = 0
    angmax = 0.0E+00

    do i = 1, nu

      if ( i /= iq .and. ( ux(i) /= xq .or. uy(i) /= yq ) ) then

        ai = angle_deg_2d ( xp, yp, xq, yq, ux(i), uy(i) )

        if ( ir == 0 .or. ai > angmax ) then

          ir = i
          xr = ux(ir)
          yr = uy(ir)
          angmax = ai
!
!  In case of ties, choose the nearer point.
!
        else if ( ir /= 0 .and. ai == angmax ) then

          di = enorm0_2d ( xq, yq, ux(i), uy(i) )
          dr = enorm0_2d ( xq, yq, xr, yr )

          if ( di < dr ) then
            ir = i
            xr = ux(ir)
            yr = uy(ir)
            angmax = ai
          end if

        end if

      end if

    end do
!
!  If we've returned to our starting point, exit.
!
    if ( ir == istart ) then
      exit
    end if

    nv = nv + 1

    if ( nv > nu ) then
      write ( *, * ) ' '
      write ( *, * ) 'POINTS_CONVEX_HULL_NLOGH_2D - Fatal error!'
      write ( *, * ) '  The algorithm failed.'
      stop
    end if
!
!  Set Q := P, P := R, and repeat.
!
    v(nv) = ir
    xp = xq
    yp = yq
    iq = ir
    xq = xr
    yq = yr

  end do
!
!  Reverse the order.
!
  call ivec_reverse ( nv, v )

  return
end
subroutine points_convex_hull_nlogn_2d ( nu, ux, uy, nv, v )
!
!*******************************************************************************
!
!! POINTS_CONVEX_HULL_NLOGN_2D computes the convex hull of 2D points.
!
!
!  Discussion:
!
!    The algorithm used requires time that is of order N * Log ( N ) in
!    N. the number of input points.
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, pages 2-8.
!
!  Modified:
!
!    30 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NU, the number of points.
!
!    Input, real UX(NU), UY(NU), the coordinates of the points.
!
!    Output, integer NV, the number of vertices in the convex hull.
!
!    Output, integer V(NU), the NV vertices that form the convex hull,
!    in counter clockwise order.   (The arrays should be dimensioned to NU.  The value
!    of NV is guaranteed to be no greater than NU.)
!
  implicit none
!
  integer nu
!
  real angle_rad_2d
  integer i
  integer l(nu+100)
  integer n1
  integer n2
  integer n3
  integer nv
  real, parameter :: pi = 3.14159265E+00
  real ux(nu)
  real uy(nu)
  integer v(nu)
  real x1
  real x2
  real x3
  real y1
  real y2
  real y3
!
!  Sort the (X,Y) points lexicographically.
!
  call rvec2_sort_a ( nu, ux, uy )
!
!  Start with L1 = ( 1, 2 )
!
!  List L1 is stored in indices N1 through N2 of L.
!
  n1 = 1
  n2 = 2

  l(1) = 1
  x2 = ux(1)
  y2 = uy(1)
  l(2) = 2
  x3 = ux(2)
  y3 = uy(2)
!
!  Append I to L1.
!    Then, while L1 contains at least 3 points
!     if the last 3 points do not make a strict right turn, delete middle one.
!
  do i = 3, nu
    n2 = n2 + 1
    l(n2) = i
    x1 = x2
    y1 = y2
    x2 = x3
    y2 = y3
    x3 = ux(i)
    y3 = uy(i)
    do while ( n2 + 1 - n1 >= 3 )

      if ( angle_rad_2d ( x1, y1, x2, y2, x3, y3 ) < pi ) then
        exit
      end if
      l(n2-1) = l(n2)
      n2 = n2 - 1
      x2 = x1
      y2 = y1
      if ( n2-2 >= 1 ) then
        x1 = ux(l(n2-2))
        y1 = uy(l(n2-2))
      end if
    end do
  end do
!
!  Start with L2 = ( N, N-1 )
!
!  List L2 is stored in indices N2 through N3 of L.
!
  n3 = n2
  n3 = n3 + 1
  x2 = ux(nu)
  y2 = uy(nu)
  l(n3) = nu - 1
  x3 = ux(nu-1)
  y3 = uy(nu-1)
!
!  Append I to L2.
!    Then, while L2 contains at least 3 points,
!    if the last 3 points do not make a strict right turn, delete middle one.
!
  do i = nu-2, 1, -1
    n3 = n3 + 1
    l(n3) = i
    x1 = x2
    y1 = y2
    x2 = x3
    y2 = y3
    x3 = ux(i)
    y3 = uy(i)
    do while ( n3 + 1 - n2 >= 3 )
      if ( angle_rad_2d ( x1, y1, x2, y2, x3, y3 ) < pi ) then
        exit
      end if
      l(n3-1) = l(n3)
      n3 = n3 - 1
      x2 = x1
      y2 = y1
      if ( n3-2 >= n2 ) then
        x1 = ux(l(n3-2))
        y1 = uy(l(n3-2))
      end if
    end do
  end do
!
!  Last entry of L is 1, so knock it off.
!
  nv = n3 - 1
  v(1:nv) = l(1:nv)
!
!  Reverse the order.
!
  call ivec_reverse ( nv, v )

  return
end
subroutine points_minidisc1_2d ( n, px, py, qx, qy, r, cx, cy )
!
!*******************************************************************************
!
!! POINTS_MINIDISC1_2D finds the smallest circle through Q containing points P.
!
!
!  Discussion:
!
!    It is assumed that there IS at least one disk which has the
!    points Q on its boundary, and which contains all the
!    points P.  For arbitrary points, this need not be the case.
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, pages 86-90.
!
!  Modified:
!
!    27 September 2000
!
!  Parameters:
!
!    Input, integer N, the number of points in the set P.
!
!    Input, real PX(N), PY(N), the X and Y coordinates of a set of
!    points in the plane.
!
!    Input, real QX, QY, a point in the plane.
!
!    Output, real R, CX, CY, the radius and center of the smallest
!    disk that encloses P and has Q on its boundary.
!
  implicit none
!
  integer n
!
  real cx
  real cy
  logical inside
  integer j
  real px(n)
  real py(n)
  real qx
  real qy
  real r
!
!  Determine the smallest disk with Q and P(1) on its boundary,
!  which is simply the circle whose diameter is the segment
!  from Q to P(1).
!
  call circle_dia2imp_2d ( qx, qy, px(1), py(1), r, cx, cy )
!
!  Now consider a point in the set P.  If it is not already
!  contained in the current circle, then expand the circle.
!
  do j = 2, n

    call circle_imp_contains_point_2d ( r, cx, cy, px(j), py(j), inside )

    if ( .not. inside ) then
      call points_minidisc2_2d ( j-1, px, py, px(j), py(j), qx, qy, r, cx, cy )
    end if

  end do

  return
end
subroutine points_minidisc2_2d ( n, px, py, q1x, q1y, q2x, q2y, r, cx, cy )
!
!*******************************************************************************
!
!! POINTS_MINIDISC2_2D finds the smallest circle through Q1 and Q2 containing points P.
!
!
!  Discussion:
!
!    It is assumed that there IS at least one disk which has the
!    points Q1 and Q2 on its boundary, and which contains all the
!    points P.  For arbitrary points, this need not be the case.
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, pages 86-90.
!
!  Modified:
!
!    27 September 2000
!
!  Parameters:
!
!    Input, integer N, the number of points in the set P.
!
!    Input, real PX(N), PY(N), the X and Y coordinates of a set of
!    points in the plane.
!
!    Input, real Q1X, Q1Y, Q2X, Q2Y, two points in the plane.
!
!    Output, real R, CX, CY, the radius and center of the smallest
!    disk that encloses P and has Q1 and Q2 on its boundary.
!
  implicit none
!
  integer n
!
  real cx
  real cy
  logical inside
  integer k
  real px(n)
  real py(n)
  real q1x
  real q1y
  real q2x
  real q2y
  real r
!
!  Determine the smallest disk with Q1, Q2 on its boundary,
!  which is simply the circle whose diameter is the segment
!  from Q1 to Q2.
!
  call circle_dia2imp_2d ( q1x, q1y, q2x, q2y, r, cx, cy )
!
!  Now consider a point in the set P.  If it is not already
!  contained in the current circle, then expand the circle.
!
  do k = 1, n

    call circle_imp_contains_point_2d ( r, cx, cy, px(k), py(k), inside )

    if ( .not. inside ) then
      call circle_exp2imp_2d ( q1x, q1y, q2x, q2y, px(k), py(k), r, cx, cy )
    end if

  end do

  return
end
subroutine points_minidisc_2d ( n, px, py, r, cx, cy )
!
!*******************************************************************************
!
!! POINTS_MINIDISC_2D finds the smallest circle containing points P.
!
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, pages 86-90.
!
!  Modified:
!
!    27 September 2000
!
!  Parameters:
!
!    Input, integer N, the number of points in the set P.
!
!    Input, real PX(N), PY(N), the X and Y coordinates of a set of
!    points in the plane.
!
!    Output, real R, CX, CY, the radius and center of the smallest
!    disk that encloses P.
!
  implicit none
!
  integer n
!
  real cx
  real cy
  integer i
  logical inside
  real px(n)
  real py(n)
  real r
!
!  N = 1
!
  if ( n == 1 ) then
    r = 0.0E+00
    cx = px(1)
    cy = py(1)
    return
  end if
!
!  N = 2
!
  if ( n == 2 ) then
    r = 0.5E+00 * sqrt ( ( px(1) - px(2) )**2 + ( py(1) - py(2) )**2 )
    cx = 0.5E+00 * ( px(1) + px(2) )
    cy = 0.5E+00 * ( py(1) + py(2) )
    return
  end if
!
!  Determine the smallest disk with P(1), P(2) on its boundary,
!
  call circle_dia2imp_2d ( px(1), py(1), px(2), py(2), r, cx, cy )
!
!  Now consider a point in the set P.  If it is not already
!  contained in the current circle, then expand the circle.
!
  do i = 3, n

    call circle_imp_contains_point_2d ( r, cx, cy, px(i), py(i), inside )

    if ( .not. inside ) then
      call points_minidisc1_2d ( i-1, px, py, px(i), py(i), r, cx, cy )
    end if

  end do

  return
end
subroutine poly_triangulate_2d ( n, x, y, triang )
!
!*******************************************************************************
!
!! POLY_TRIANGULATE_2D returns a triangulation of a polygon.
!
!
!  Discussion:
!
!    Given a polygon with N (distinct) nodes, a triangulation consists of
!    N-2 (distinct) triangles, with the property that each triangle is
!    constructed from vertices of the polygon, and that every point in the
!    polygon is contained in exactly one triangle, unless it lies on the
!    boundary of two or more triangles, and that no point exterior to the
!    polygon is contained in any triangle.  In other words, a triangulation
!    is a dissection of the area of the polygon into triangles.
!
!    This algorithm should be correct, but not efficient.
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, pages 46-49.
!
!  Modified:
!
!    09 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of nodes in the polygon.
!
!    Input, real X(N), Y(N), the coordinates of the nodes, listed in
!    counter-clockwise order.
!
!    Output, integer TRIANG(3,N-2), the triangulation of the polygon.
!
  implicit none
!
  integer n
  integer, parameter :: stack2_max = 100
!
  integer best
  integer degree
  integer degree2
  real dist
  real dist_max
  integer i
  logical inside
  integer ipol
  integer j
  integer number
  integer poly(n)
  integer stack1(n-2)
  integer stack1_num
  integer stack2(stack2_max)
  integer stack2_num
  integer t
  integer triang(3,n-2)
  integer triang_num
  integer u
  integer v
  integer w
  real x(n)
  real x1
  real x2
  real x3
  real x4
  real y(n)
  real y1
  real y2
  real y3
  real y4
!
  if ( n <= 2 ) then
    return
  end if

  degree = n
  call ivec_identity ( n, poly )

  call poly_reorder_nodes ( n, x, y, degree, poly )

  stack1_num = 0
  stack2_num = 0
  triang_num = 0

  do

    if ( degree == 3 ) then

      triang_num = triang_num + 1
      triang(1:3,triang_num) = poly(1:3)

      if ( stack1_num <= 0 ) then
        exit
      end if

      call ivec_pop ( degree, poly, n-2, stack1_num, stack1, &
        stack2_max, stack2_num, stack2 )

      call poly_reorder_nodes ( n, x, y, degree, poly )

      cycle

    end if

    u = poly(degree)
    v = poly(1)
    w = poly(2)

    x1 = x(u)
    y1 = y(u)
    x2 = x(v)
    y2 = y(v)
    x3 = x(w)
    y3 = y(w)

    dist_max = 0.0E+00
    best = 0
    number = 0

    do j = 3, degree-1

      t = j
      x4 = x(poly(t))
      y4 = y(poly(t))

      call triangle_contains_point_2d ( x1, y1, x2, y2, x3, y3, x4, y4, inside )

      if ( inside ) then

        number = number + 1
        call line_exp_point_dist_2d ( x1, y1, x3, y3, x4, y4, dist )

        if ( dist >= dist_max ) then
          dist_max = dist
          best = t
        end if

      end if

    end do
!
!  If there were no polygonal nodes inside the triangle (U,V,W),
!  then triangle (U,V,W) can be used as part of the triangulation,
!  and the remaining region is simply the polygon (W,...,U).
!
    if ( number == 0 ) then

      triang_num = triang_num + 1
      triang(1,triang_num) = poly(1)
      triang(2,triang_num) = poly(2)
      triang(3,triang_num) = poly(degree)

      poly(1:n-1) = poly(2:n)
      degree = degree - 1

      call poly_reorder_nodes ( n, x, y, degree, poly )

    else
!
!  If there were polygonal nodes inside the triangle (U,V,W),
!  then POLY(T) is the node that was furthest from the line (U,W).
!
!  Split the polygon (V,W,part1,T,part2,U)
!  into (V,W,part1,T) and (V,T,part2,U)
!
      degree2 = t

      call ivec_push ( degree2, poly(1:t), n-2, stack1_num, stack1, &
        stack2_max, stack2_num, stack2 )

      degree2 = degree + 2 - t
      poly(2:degree2) = poly(t:degree)
      degree = degree2

      call poly_reorder_nodes ( n, x, y, degree, poly )

    end if

  end do

  return
end
subroutine poly_reorder_nodes ( nxy, x, y, npoly, poly )
!
!*******************************************************************************
!
!! POLY_REORDER_NODES reorders nodes of a polygon so node 1 is leftest lowest.
!
!
!  Modified:
!
!    12 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NXY, the number of nodes.
!
!    Input, real X(NXY), Y(NXY), the coordinates of the nodes.
!
!    Input, integer NPOLY, the number of nodes of the polygon.
!
!    Input/output, POLY(NPOLY), the indices of the nodes.
!
  implicit none
!
  integer npoly
  integer nxy
!
  integer i
  integer imin
  integer p
  integer pmin
  integer poly(npoly)
  integer poly2(npoly)
  real x(nxy)
  real y(nxy)
!
  imin = 1
  pmin = poly(imin)

  do i = 2, npoly
    p = poly(i)
    if ( &
      ( x(p) < x(pmin) ) .or. &
      ( x(p) == x(pmin) .and. y(p) < y(pmin) ) ) then
      imin = i
      pmin = p
    end if
  end do

  if ( imin /= 1 ) then

    poly2(1:npoly+1-imin) = poly(imin:npoly)
    poly2(npoly+2-imin:npoly) = poly(1:imin-1)

    poly(1:npoly) = poly2(1:npoly)

  end if

  return
end
subroutine polycon_minkowski_sum_linear ( nu, ux, uy, nv, vx, vy, nw, wx, wy )
!
!*******************************************************************************
!
!! POLYCON_MINKOWSKI_SUM_LINEAR computes the Minkowski sum of two convex polygons.
!
!
!  Discussion:
!
!    For two geometric shapes U and V, the Minkowski sum W is the
!    set of all points
!      w = u + v
!    formed by adding points from the two shapes.
!
!    The Minkowski sum of two convex polygons is also a convex polygon.
!
!    The algorithm used here is only valid for convex polygons.
!
!    The vertices must be listed in counterclockwise order, with
!    the first vertex having the smallest Y coordinate.  If two
!    or more vertices have the same Y coordinate, then the one with
!    smallest X coordinate should be first.
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, pages 275-281.
!
!  Modified:
!
!    30 August 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NU, the number of vertices of the first polygon.
!
!    Input, real UX(NU), UY(NU), the coordinates of the vertices of the
!    first polygon.
!
!    Input, integer NV, the number of vertices of the second polygon.
!
!    Input, real VX(NV), VY(NY), the coordinates of the vertices of the
!    second polygon.
!
!    Output, integer NW, the number of vertices of the sum polygon.
!    NW will be at most NV+NU.
!
!    Output, real WX(*), WY(*), the coordinates of the vertices of the
!    sum polygon.
!
  implicit none
!
  integer nu
  integer nv
!
  real angle_rad_2d
  integer i
  integer i_wrap
  integer iw
  integer ip1w
  integer j
  integer jp1w
  integer jw
  integer nw
  real pi
  real ux(nu)
  real uy(nu)
  real u_angle
  logical u_done
  real vx(nv)
  real vy(nv)
  real v_angle
  logical v_done
  real wx(nu+nv+2)
  real wy(nu+nv+2)
!
  i = 1
  iw = 1
  u_done = .false.
  j = 1
  jw = 1
  v_done = .false.
  nw = 0

  do

    nw = nw + 1

    wx(nw) = ux(iw) + vx(jw)
    wy(nw) = uy(iw) + vy(jw)

    ip1w = i_wrap ( i+1, 1, nu )
    u_angle = angle_rad_2d ( ux(iw) + 1.0E+00, uy(iw), ux(iw), uy(iw), &
      ux(ip1w), uy(ip1w) )
    u_angle = u_angle + real ( ( i - nu ) / nu ) * 2.0E+00 * pi ()

    jp1w = i_wrap ( j+1, 1, nv )
    v_angle = angle_rad_2d ( vx(jw) + 1.0E+00, vy(jw), vx(jw), vy(jw), &
      vx(jp1w), vy(jp1w) )
    v_angle = v_angle + real ( ( j - nv ) / nv ) * 2.0E+00 * pi ()

    if ( u_angle < v_angle ) then
      i = i + 1
    else if ( u_angle > v_angle ) then
      j = j + 1
    else
      i = i + 1
      j = j + 1
    end if

    if ( i == nu + 1 ) then
      u_done = .true.
    end if

    if ( j == nv + 1 ) then
      v_done = .true.
    end if

    if ( u_done .and. v_done ) then
      exit
    end if

    iw = i_wrap ( i, 1, nu )
    jw = i_wrap ( j, 1, nv )

  end do

  return
end
subroutine polycon_minkowski_sum_n2logn2 ( nu, ux, uy, nv, vx, vy, nw, wx, &
  wy )
!
!*******************************************************************************
!
!! POLYCON_MINKOWSKI_SUM_N2LOGN2 Minkowski sums two convex polygons.
!
!
!  Discussion:
!
!    For two geometric shapes U and V, the Minkowski sum W is the
!    set of all points
!      w = u + v
!    formed by adding points from the two shapes.
!
!    The Minkowski sum of two convex polygons is also a convex polygon.
!
!    The algorithm used here is only valid for convex polygons.
!
!    The vertices must be listed in counterclockwise order, with
!    the first vertex having the smallest Y coordinate.  If two
!    or more vertices have the same Y coordinate, then the one with
!    smallest X coordinate should be first.
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, pages 275-281.
!
!  Modified:
!
!    11 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NU, the number of vertices of the first polygon.
!
!    Input, real UX(NU), UY(NU), the coordinates of the vertices of the
!    first polygon.
!
!    Input, integer NV, the number of vertices of the second polygon.
!
!    Input, real VX(NV), VY(NY), the coordinates of the vertices of the
!    second polygon.
!
!    Output, integer NW, the number of vertices of the sum polygon.
!    NW will be at most NU+NV.
!
!    Output, real WX(*), WY(*), the coordinates of the vertices of the
!    sum polygon.
!
  implicit none
!
  integer nu
  integer nv
!
  integer i
  integer j
  integer nuv
  integer nw
  real uvx(nu*nv)
  real uvy(nu*nv)
  real ux(nu)
  real uy(nu)
  real vx(nv)
  real vy(nv)
  integer w(nu+nv)
  real wx(nu+nv)
  real wy(nu+nv)
!
!  Generate points from all pairs.
!
  nuv = 0
  do i = 1, nu
    do j = 1, nv
      nuv = nuv + 1
      uvx(nuv) = ux(i) + vx(j)
      uvy(nuv) = uy(i) + vy(j)
    end do
  end do
!
!  Compute the convex hull.
!  The output value of NW should be no more than NU+NV.
!
  call points_convex_hull_nlogn_2d ( nuv, uvx, uvy, nw, w )
!
!  Collect the points.
!
  do i = 1, nw
    wx(i) = uvx(w(i))
    wy(i) = uvy(w(i))
  end do

  return
end
subroutine r_swap ( x, y )
!
!*******************************************************************************
!
!! R_SWAP swaps two real values.
!
!
!  Modified:
!
!    01 May 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real X, Y.  On output, the values of X and
!    Y have been interchanged.
!
  implicit none
!
  real x
  real y
  real z
!
  z = x
  x = y
  y = z

  return
end
function radians_to_degrees ( angle )
!
!*******************************************************************************
!
!! RADIANS_TO_DEGREES converts an angle from radians to degrees.
!
!
!  Modified:
!
!    10 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ANGLE, an angle in radians.
!
!    Output, real RADIANS_TO_DEGREES, the equivalent angle
!    in degrees.
!
  implicit none
!
  real angle
  real, parameter :: pi = &
    3.14159265358979323846264338327950288419716939937510E+00
  real radians_to_degrees
!
  radians_to_degrees = ( angle / pi ) * 180.0E+00

  return
end
subroutine rect_int_2d ( x1, y1, x2, y2, x3, y3, x4, y4, flag, x5, y5, x6, y6 )
!
!*******************************************************************************
!
!! RECT_INT_2D computes the intersection of two rectangles in 2D.
!
!
!  Modified:
!
!    07 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, the corners of the first segment.
!
!    Input, real X3, Y3, X4, Y4, the corners of the second segment.
!
!    Output, integer FLAG, records the results.
!    0, the rectangles do not intersect.
!    1, the intersection is a point.
!    2, the intersection is a line.
!    3, the intersection is a rectangle.
!
!    Output, real X5, Y5, X6, Y6, the corners of the intersection.
!    If FLAG = 0, X5 = Y5 = X6 = Y6 = 0.
!
  implicit none
!
  integer flag
  real x1
  real x2
  real x3
  real x4
  real x5
  real x6
  real y1
  real y2
  real y3
  real y4
  real y5
  real y6
!
  x5 = 0.0E+00
  y5 = 0.0E+00
  x6 = 0.0E+00
  y6 = 0.0E+00

  call lines_seg_int_1d ( x1, x2, x3, x4, flag, x5, x6 )

  if ( flag == 0 ) then
    return
  end if

  call lines_seg_int_1d ( y1, y2, y3, y4, flag, y5, y6 )

  return
end
subroutine rmat_print ( lda, m, n, a, title )
!
!*******************************************************************************
!
!! RMAT_PRINT prints a real matrix.
!
!
!  Modified:
!
!    24 March 2000
!
!  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, real A(LDA,N), the matrix to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none
!
  integer lda
  integer n
!
  real a(lda,n)
  integer i
  integer j
  integer jhi
  integer jlo
  integer m
  character ( len = * ) title
!
  if ( title /= ' ' ) then
    write ( *, * ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  do jlo = 1, n, 5
    jhi = min ( jlo + 4, n )
    write ( *, * ) ' '
    write ( *, '(6x,5(i7,7x))' ) ( j, j = jlo, jhi )
    write ( *, * ) ' '
    do i = 1, m
      write ( *, '(i6,5g14.6)' ) i, a(i,jlo:jhi)
    end do
  end do

  return
end
subroutine rmat_solve ( a, n, nrhs, info )
!
!*******************************************************************************
!
!! RMAT_SOLVE uses Gauss-Jordan elimination to solve an N by N linear system.
!
!
!  Modified:
!
!    08 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real A(N,N+NRHS), contains in rows and columns 1
!    to N the coefficient matrix, and in columns N+1 through
!    N+NRHS, the right hand sides.  On output, the coefficient matrix
!    area has been destroyed, while the right hand sides have
!    been overwritten with the corresponding solutions.
!
!    Input, integer NRHS, the number of right hand sides.  NRHS
!    must be at least 0.
!
!    Output, integer INFO, singularity flag.
!    0, the matrix was not singular, the solutions were computed;
!    J, factorization failed on step J, and the solutions could not
!    be computed.
!
  implicit none
!
  integer n
  integer nrhs
!
  real a(n,n+nrhs)
  real apivot
  real factor
  integer i
  integer info
  integer ipivot
  integer j
  integer k
  real temp
!
  info = 0

  do j = 1, n
!
!  Choose a pivot row.
!
    ipivot = j
    apivot = a(j,j)

    do i = j+1, n
      if ( abs ( a(i,j) ) > abs ( apivot ) ) then
        apivot = a(i,j)
        ipivot = i
      end if
    end do

    if ( apivot == 0.0E+00 ) then
      info = j
      return
    end if
!
!  Interchange.
!
    do i = 1, n + nrhs
      call r_swap ( a(ipivot,i), a(j,i) )
    end do
!
!  A(J,J) becomes 1.
!
    a(j,j) = 1.0E+00
    a(j,j+1:n+nrhs) = a(j,j+1:n+nrhs) / apivot
!
!  A(I,J) becomes 0.
!
    do i = 1, n

      if ( i /= j ) then

        factor = a(i,j)
        a(i,j) = 0.0E+00
        a(i,j+1:n+nrhs) = a(i,j+1:n+nrhs) - factor * a(j,j+1:n+nrhs)

      end if

    end do

  end do

  return
end
subroutine rmat2_inverse ( a, b, det )
!
!*******************************************************************************
!
!! RMAT2_INVERSE inverts a 2 by 2 real matrix using Cramer's rule.
!
!
!  Modified:
!
!    16 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real A(2,2), the matrix to be inverted.
!
!    Output, real B(2,2), the inverse of the matrix A.
!
!    Output, real DET, the determinant of the matrix A.
!
!    If DET is zero, then A is singular, and does not have an
!    inverse.  In that case, B is simply set to zero, and a
!    message is printed.
!
!    If DET is nonzero, then its value is roughly an estimate
!    of how nonsingular the matrix A is.
!
  implicit none
!
  real a(2,2)
  real b(2,2)
  real det
!
!  Compute the determinant.
!
  det = a(1,1) * a(2,2) - a(1,2) * a(2,1)
!
!  If the determinant is zero, bail out.
!
  if ( det == 0.0E+00 ) then

    b(1:2,1:2) = 0.0E+00

    return
  end if
!
!  Compute the entries of the inverse matrix using an explicit formula.
!
  b(1,1) = + a(2,2) / det
  b(1,2) = - a(1,2) / det
  b(2,1) = - a(2,1) / det
  b(2,2) = + a(1,1) / det

  return
end
subroutine rvec2_compare ( n, a1, a2, i, j, isgn )
!
!*******************************************************************************
!
!! RVEC2_COMPARE compares pairs of reals stored in two vectors.
!
!
!  Discussion:
!
!    The lexicographic ordering is used.
!
!  Modified:
!
!    08 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data items.
!
!    Input, real A1(N), A2(N), contain the two components of each item.
!
!    Input, integer I, J, the items to be compared.
!
!    Output, integer ISGN, the results of the comparison:
!    -1, item I < item J,
!     0, item I = item J,
!    +1, item I > item J.
!
  implicit none
!
  integer n
!
  real a1(n)
  real a2(n)
  integer i
  integer isgn
  integer j
!
  isgn = 0

       if ( a1(i) < a1(j) ) then

    isgn = -1

  else if ( a1(i) == a1(j) ) then

         if ( a2(i) < a2(j) ) then
      isgn = -1
    else if ( a2(i) < a2(j) ) then
      isgn = 0
    else if ( a2(i) > a2(j) ) then
      isgn = +1
    end if

  else if ( a1(i) > a1(j) ) then

    isgn = +1

  end if

  return
end
subroutine rvec2_print ( n, a1, a2, title )
!
!*******************************************************************************
!
!! RVEC2_PRINT prints a pair of real vectors.
!
!
!  Modified:
!
!    14 June 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, real A1(N), A2(N), the vectors to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none
!
  integer n
!
  real a1(n)
  real a2(n)
  integer i
  character ( len = * ) title
!
  if ( title /= ' ' ) then
    write ( *, * ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, * ) ' '
  do i = 1, n
    write ( *, '(i6,2g14.6)' ) i, a1(i), a2(i)
  end do

  return
end
subroutine rvec2_sort_a ( n, a1, a2 )
!
!*******************************************************************************
!
!! RVEC2_SORT_A ascending sorts a vector of pairs of integers.
!
!
!  Discussion:
!
!    Each item to be sorted is a pair (I,J), with the I
!    and J values stored in separate vectors A1 and A2.
!
!  Modified:
!
!    08 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of items of data.
!
!    Input/output, real A1(N), A2(N), the data to be sorted.
!
  implicit none
!
  integer n
!
  real a1(n)
  real a2(n)
  integer i
  integer indx
  integer isgn
  integer j
!
!  Initialize.
!
  i = 0
  indx = 0
  isgn = 0
  j = 0
!
!  Call the external heap sorter.
!
  do

    call sort_heap_external ( n, indx, i, j, isgn )
!
!  Interchange the I and J objects.
!
    if ( indx > 0 ) then

      call r_swap ( a1(i), a1(j) )
      call r_swap ( a2(i), a2(j) )
!
!  Compare the I and J objects.
!
    else if ( indx < 0 ) then

      call rvec2_compare ( n, a1, a2, i, j, isgn )

    else if ( indx == 0 ) then

      exit

    end if

  end do

  return
end
subroutine rvec_reverse ( n, a )
!
!*******************************************************************************
!
!! RVEC_REVERSE reverses the elements of a real vector.
!
!
!  Example:
!
!    Input:
!
!      N = 5, A = ( 11.0, 12.0, 13.0, 14.0, 15.0E+00 ).
!
!    Output:
!
!      A = ( 15.0, 14.0, 13.0, 12.0, 11.0E+00 ).
!
!  Modified:
!
!    06 October 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the array.
!
!    Input/output, real A(N), the array to be reversed.
!
  implicit none
!
  integer n
!
  real a(n)
  integer i
!
  do i = 1, n/2
    call r_swap ( a(i), a(n+1-i) )
  end do

  return
end
subroutine sort_heap_external ( n, indx, i, j, isgn )
!
!*******************************************************************************
!
!! SORT_HEAP_EXTERNAL externally sorts a list of items into linear order.
!
!
!  Discussion:
!
!    The actual list of data is not passed to the routine.  Hence this
!    routine may be used to sort integers, reals, numbers, names,
!    dates, shoe sizes, and so on.  After each call, the routine asks
!    the user to compare or interchange two items, until a special
!    return value signals that the sorting is completed.
!
!  Modified:
!
!    12 November 2000
!
!  Reference:
!
!    A Nijenhuis and H Wilf,
!    Combinatorial Algorithms,
!    Academic Press, 1978, second edition,
!    ISBN 0-12-519260-6.
!
!  Parameters:
!
!    Input, integer N, the number of items to be sorted.
!
!    Input/output, integer INDX, the main communication signal.
!
!    The user must set INDX to 0 before the first call.
!    Thereafter, the user should not change the value of INDX until
!    the sorting is done.
!
!    On return, if INDX is
!
!      greater than 0,
!      * interchange items I and J;
!      * call again.
!
!      less than 0,
!      * compare items I and J;
!      * set ISGN = -1 if I precedes J, ISGN = +1 if J precedes I;
!      * call again.
!
!      equal to 0, the sorting is done.
!
!    Output, integer I, J, the indices of two items.
!    On return with INDX positive, elements I and J should be interchanged.
!    On return with INDX negative, elements I and J should be compared, and
!    the result reported in ISGN on the next call.
!
!    Input, integer ISGN, results of comparison of elements I and J.
!    (Used only when the previous call returned INDX less than 0).
!    ISGN <= 0 means I precedes J;
!    ISGN => 0 means J precedes I.
!
  implicit none
!
  integer i
  integer indx
  integer isgn
  integer j
  integer, save :: k = 0
  integer, save :: k1 = 0
  integer n
  integer, save :: n1 = 0
!
!  INDX = 0: This is the first call.
!
  if ( indx == 0 ) then

    n1 = n
    k = n / 2
    k1 = k
!
!  INDX < 0: The user is returning the results of a comparison.
!
  else if ( indx < 0 ) then

    if ( indx == -2 ) then

      if ( isgn < 0 ) then
        i = i + 1
      end if

      j = k1
      k1 = i
      indx = - 1
      return

    end if

    if ( isgn > 0 ) then
      indx = 2
      return
    end if

    if ( k <= 1 ) then

      if ( n1 == 1 ) then
        indx = 0
      else
        i = n1
        n1 = n1 - 1
        j = 1
        indx = 1
      end if

      return

    end if

    k = k - 1
    k1 = k
!
!  INDX > 0, the user was asked to make an interchange.
!
  else if ( indx == 1 ) then

    k1 = k

  end if

  do

    i = 2 * k1

    if ( i == n1 ) then
      j = k1
      k1 = i
      indx = - 1
      return
    else if ( i <= n1 ) then
      j = i + 1
      indx = - 2
      return
    end if

    if ( k <= 1 ) then
      exit
    end if

    k = k - 1
    k1 = k

  end do

  if ( n1 == 1 ) then
    indx = 0
  else
    i = n1
    n1 = n1 - 1
    j = 1
    indx = 1
  end if

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

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

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

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

  return
end
subroutine triangle_contains_point_2d ( x1, y1, x2, y2, x3, y3, x, y, inside )
!
!*******************************************************************************
!
!! TRIANGLE_CONTAINS_POINT_2D finds if a point is inside a triangle in 2D.
!
!
!  Modified:
!
!    12 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real X1, Y1, X2, Y2, X3, Y3, the coordinates of
!    the corners of the triangle.
!
!    Input, real X, Y, the point to be checked.
!
!    Output, logical INSIDE, is .TRUE. if (X,Y) is inside
!    the triangle or on its boundary, and .FALSE. otherwise.
!
  implicit none
!
  integer, parameter :: N = 2
  integer, parameter :: NRHS = 1
!
  real a(N,N+NRHS)
  real c1
  real c2
  integer info
  logical inside
  real x
  real x1
  real x2
  real x3
  real y
  real y1
  real y2
  real y3
!
!  Set up the linear system
!
!    ( X2-X1  X3-X1 ) C1  = X-X1
!    ( Y2-Y1  Y3-Y1 ) C2    Y-Y1
!
!  which is satisfied by the barycentric coordinates of (X,Y).
!
  a(1,1) = x2 - x1
  a(1,2) = x3 - x1
  a(1,3) = x - x1

  a(2,1) = y2 - y1
  a(2,2) = y3 - y1
  a(2,3) = y - y1
!
!  Solve the linear system.
!
  call rmat_solve ( a, N, NRHS, info )

  if ( info /= 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'TRIANGLE_CONTAINS_POINT - Fatal error!'
    write ( *, * ) '  The linear system is singular.'
    write ( *, * ) '  The input data does not form a proper triangle.'
    stop
  end if

  c1 = a(1,3)
  c2 = a(2,3)
!
!  If the point is in the triangle, its barycentric coordinates
!  must both be nonnegative, and sum to no more than 1.
!
  if ( c1 < 0.0E+00 .or. c2 < 0.0E+00 ) then
    inside = .false.
  else if ( c1 + c2 > 1.0E+00 ) then
    inside = .false.
  else
    inside = .true.
  end if

  return
end
subroutine triangulate_tricolor ( node_num, triang, color )
!
!*******************************************************************************
!
!! TRIANGULATE_TRICOLOR three-colors the nodes of a triangulated polygon.
!
!
!  Discussion:
!
!    While the data in TRIANG must represent a proper triangulation
!    of the polygon, no check for this is made.
!
!    The three-coloring is done in such a way that every triangle
!    of the triangulation has one node of each color.
!
!    This solves the art gallery problem, which asks for the minimum
!    number of guards necessary in a (polygonal) art gallery, so that
!    every point inside the gallery is seen by at least one guard.
!    The answer is that (N/3) (integer division, with truncation) guards
!    are sufficient.  Triangulate the polygon, tricolor the triangulation,
!    choose the color that shows up least often (in case N is not divisible
!    by 3) and place a guard at those points.  Since every triangle has
!    a node of that color, and the triangles make up the gallery, you're done.
!    This assumes that the gallery is a simple polygon.
!
!  Reference:
!
!    de Berg, van Krevald, Overmars, Schwarzkopf,
!    Computational Geometry,
!    Second Edition,
!    Springer, 2000, pages 47-48.
!
!  Modified:
!
!    17 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NODE_NUM, the number of nodes in the polygon.
!
!    Input, integer TRIANG(3,NODE_NUM-2), the triangulation of the polygon.
!
!    Output, integer COLOR(NODE_NUM), an assignment of the "colors" 1, 2
!    and 3 to the triangles in such a way that every triangle in the
!    triangulation has one node of each color.
!
  implicit none
!
  integer node_num
!
  integer color(node_num)
  integer color1
  integer color2
  integer color3
  integer node1
  integer node2
  integer node3
  integer stack(6*(node_num-3))
  integer stack_max
  integer stack_num
  integer t1
  integer t2
  integer triang(3,node_num-2)
!
  stack_max = 6 * ( node_num - 3 )
  t1 = 1

  node1 = triang(1,t1)
  node2 = triang(2,t1)
  node3 = triang(3,t1)

  color1 = 1
  color2 = 2
  color3 = 3

  color(node1) = color1
  color(node2) = color2
  color(node3) = color3

  stack_num = 0

  call triangulate_color_push ( t1, node1, color1, node2, color2, color3, &
    stack_max, stack_num, stack )

  call triangulate_color_push ( t1, node2, color2, node3, color3, color1, &
    stack_max, stack_num, stack )

  call triangulate_color_push ( t1, node3, color3, node1, color1, color2, &
    stack_max, stack_num, stack )

  do while ( stack_num > 0 )

    call triangulate_color_pop ( t1, node1, color1, node2, color2, color3, &
      stack_max, stack_num, stack )

    call triangulate_common_edge ( node_num-2, triang, node1, node2, t1, &
      t2, node3 )

    if ( t2 /= 0 ) then

      color(node3) = color3

      call triangulate_color_push ( t2, node1, color1, node3, color3, color2, &
        stack_max, stack_num, stack )

      call triangulate_color_push ( t2, node3, color3, node2, color2, color1, &
        stack_max, stack_num, stack )

    end if

  end do

  return
end
subroutine triangulate_color_push ( t, node1, color1, node2, color2, color3, &
  stack_max, stack_num, stack )
!
!*******************************************************************************
!
!! TRIANGULATE_COLOR_PUSH pushes a side of a colored triangle onto the stack.
!
!
!  Discussion:
!
!    This is a utility routine used by TRIANGULATE_TRICOLOR.
!
!  Modified:
!
!    17 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer T, the triangle that the side belongs to.
!
!    Input, integer NODE1, COLOR1, the starting node of the edge, and its color.
!
!    Input, integer NODE2, COLOR2, the end node of the edge, and its color.
!
!    Input, integer COLOR3, the remaining color.
!
!    Input, integer STACK_MAX, the maximum size of the stack.
!
!    Input/output, integer STACK_NUM, the current size of the stack.
!
!    Input/output, integer STACK(STACK_MAX), the stack.
!
  implicit none
!
  integer stack_max
!
  integer color1
  integer color2
  integer color3
  integer node1
  integer node2
  integer stack(stack_max)
  integer stack_num
  integer t
!
  stack_num = stack_num + 1
  stack(stack_num) = t
  stack_num = stack_num + 1
  stack(stack_num) = node1
  stack_num = stack_num + 1
  stack(stack_num) = color1
  stack_num = stack_num + 1
  stack(stack_num) = node2
  stack_num = stack_num + 1
  stack(stack_num) = color2
  stack_num = stack_num + 1
  stack(stack_num) = color3

  return
end
subroutine triangulate_color_pop ( t, node1, color1, node2, color2, color3, &
  stack_max, stack_num, stack )
!
!*******************************************************************************
!
!! TRIANGULATE_COLOR_POP pops a side of a colored triangle from the stack.
!
!
!  Discussion:
!
!    This is a utility routine used by TRIANGULATE_TRICOLOR.
!
!  Modified:
!
!    17 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer T, the triangle that the side belongs to.
!
!    Output, integer NODE1, COLOR1, the starting node of the edge, and 
!    its color.
!
!    Output, integer NODE2, COLOR2, the end node of the edge, and its color.
!
!    Output, integer COLOR3, the remaining color.
!
!    Input, integer STACK_MAX, the maximum size of the stack.
!
!    Input/output, integer STACK_NUM, the current size of the stack.
!
!    Input/output, integer STACK(STACK_MAX), the stack.
!
  implicit none
!
  integer stack_max
!
  integer color1
  integer color2
  integer color3
  integer node1
  integer node2
  integer stack(stack_max)
  integer stack_num
  integer t
!
  color3 = stack(stack_num)
  stack_num = stack_num - 1
  color2 = stack(stack_num)
  stack_num = stack_num - 1
  node2 = stack(stack_num)
  stack_num = stack_num - 1
  color1 = stack(stack_num)
  stack_num = stack_num - 1
  node1 = stack(stack_num)
  stack_num = stack_num - 1
  t = stack(stack_num)
  stack_num = stack_num - 1

  return
end
subroutine triangulate_common_edge ( triang_num, triang, node1, node2, t1, &
  t2, node3 )
!
!*******************************************************************************
!
!! TRIANGULATE_COMMON_EDGE seeks the other triangle that shares an edge.
!
!
!  Discussion:
!
!    This is a utility routine used by TRIANGULATE_TRICOLOR.
!
!  Modified:
!
!    17 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TRIANG_NUM, the number of triangles.
!
!    Input, integer TRIANG(3,TRIANG_NUM), the triangulation of the polygon.
!
!    Input, integer NODE1, NODE2, the starting and ending nodes of an edge.
!
!    Input, integer T1, the triangle to which the edge (NODE1,NODE2) belongs.
!
!    Output, integer T2, is 0 if there is no triangle containing the
!    matching edge (NODE2,NODE1), or it is the index of a triangle containing
!    the edge (NODE2,NODE1).
!
!    Output, integer NODE3, the other node in triangle T2.
!
  implicit none
!
  integer triang_num
!
  integer i_wrap
  integer j1
  integer j2
  integer j3
  integer node1
  integer node2
  integer node3
  integer t
  integer t1
  integer t2
  integer triang(3,triang_num)
!
  t2 = 0
  node3 = 0

  do t = 1, triang_num
    do j1 = 1, 3
      j2 = i_wrap ( j1+1, 1, 3 )
      if ( triang(j2,t) == node1 .and. triang(j1,t) == node2 ) then
        t2 = t
        j3 = i_wrap ( j1+2, 1, 3 )
        node3 = triang(j3,t)
        return
      end if
    end do
  end do

  return
end
subroutine triangulation_boundary_count ( point_num, tri_num, bound_num )
!
!*******************************************************************************
!
!! TRIANGULATION_BOUNDARY_COUNT returns the number of boundary edges.
!
!
!  Discussion:
!
!    We assume we are given information about a legal, maximal triangulation
!    of a set of points in the plane.  In a maximal triangulation, no more
!    edges can be added without crossing an existing edge.
!
!    Given the number of points and triangles, we are going to apply
!    Euler's formula to determine the number of edges that lie on the
!    convex hull of the set of points.
!
!    The number of faces, including the infinite face, is TRI_NUM + 1.
!
!    Let BOUND_NUM denote the number of edges on the boundary of the convex
!    hull.  Each of the TRI_NUM triangles uses three edges.  Every edge
!    occurs in two different faces, so the number of edges must be
!    ( 3 * TRI_NUM + BOUND_NUM ) / 2.
!
!    The number of points is POINT_NUM.
!
!    Euler's formula asserts that, for a simple connected figure in the
!    plane with no edge crossings, POINT_NUM points, EDGE_NUM edges and
!    FACE_NUM faces:
!
!      POINT_NUM - EDGE_NUM + FACE_NUM = 2
!
!    In our context, this becomes
!
!      POINT_NUM - ( 3 * TRI_NUM + BOUND_NUM ) / 2 + TRI_NUM + 1 = 2
!
!    or
!
!      BOUND_NUM = 2 * POINT_NUM - TRI_NUM - 2
!
!  Modified:
!
!    25 July 2001
!
!  Reference:
!
!    de Berg, Krevald, Overmars and Schwarzkopf,
!    Computational Geometry, Section 9.1,
!    Springer, 2000.
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer TRI_NUM, the number of triangles.
!
!    Output, integer BOUND_NUM, the number of edges that lie on the convex
!    hull of the triangulation.
!
  implicit none
!
  integer bound_num
  integer point_num
  integer tri_num
!
  bound_num = 2 * point_num - tri_num - 2

  return
end