function angle_contains_point_2d ( p1, p2, p3, p )

!*****************************************************************************80
!
!! ANGLE_CONTAINS_POINT_2D determines if an angle contains a point, in 2D.
!
!  Discussion:
!
!    The angle is defined by the sequence of points P1, P2 and P3.
!
!    The point is "contained" by the angle if the ray P - P2
!    is between (in a counterclockwise sense) the rays P1 - P2
!    and P3 - P2.
!
!        P1
!        /
!       /   P
!      /  .
!     / .
!    P2--------->P3
!
!  Modified:
!
!    14 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) P1(2), P2(2), P3(2), three points that define 
!    the angle.  The order of these points matters!
!
!    Input, real ( kind = 8 ) P(2), the point to be checked.
!
!    Output, logical ANGLE_CONTAINS_POINT_2D, is TRUE if the point 
!    is inside the angle.
!
  implicit none

  integer, parameter :: dim_num = 2

  logical angle_contains_point_2d
  real ( kind = 8 ) angle_rad_2d
  real ( kind = 8 ) p(dim_num)
  real ( kind = 8 ) p1(dim_num)
  real ( kind = 8 ) p2(dim_num)
  real ( kind = 8 ) p3(dim_num)

  if ( angle_rad_2d ( p1, p2, p ) <= angle_rad_2d ( p1, p2, p3 ) ) then
    angle_contains_point_2d = .true.
  else
    angle_contains_point_2d = .false.
  end if

  return
end
subroutine angle_half_2d ( p1, p2, p3, p4 )

!*****************************************************************************80
!
!! ANGLE_HALF_2D finds half an angle in 2D.
!
!  Discussion:
!
!    The original angle is defined by the sequence of points P1, P2 and P3.
!
!    The point P4 is calculated so that:
!
!      (P1,P2,P4) = (P1,P2,P3) / 2
!
!        P1
!        /
!       /   P4
!      /  .
!     / .
!    P2--------->P3
!
!  Modified:
!
!    01 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) P1(2), P2(2), P3(2), points defining the angle.
!
!    Input, real ( kind = 8 ) P4(2), a point defining the half angle.
!    The vector P4 - P2 will have unit norm.
!
  implicit none

  integer, parameter :: dim_num = 2

  real ( kind = 8 ) p1(dim_num)
  real ( kind = 8 ) p2(dim_num)
  real ( kind = 8 ) p3(dim_num)
  real ( kind = 8 ) p4(dim_num)

  p4(1:2) = 0.5D+00 * ( &
      ( p1(1:2) - p2(1:2) ) / sqrt ( sum ( ( p1(1:2) - p2(1:2) )**2 ) ) &
    + ( p3(1:2) - p2(1:2) ) / sqrt ( sum ( ( p3(1:2) - p2(1:2) )**2 ) ) )

  p4(1:2) = p2(1:2) + p4(1:2) / sqrt ( sum ( p4(1:2)**2 ) )

  return
end
function angle_rad_2d ( p1, p2, p3 )

!*****************************************************************************80
!
!! 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 ( P1, P2, P3 ) + ANGLE_RAD_2D ( P3, P2, P1 ) = 2 * PI
!
!        P1
!        /
!       /
!      /
!     /
!    P2--------->P3
!
!  Modified:
!
!    14 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) P1(2), P2(2), P3(2), define the rays
!    P1 - P2 and P3 - P2 which define the angle.
!
!    Output, real ( kind = 8 ) 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

  integer, parameter :: dim_num = 2

  real ( kind = 8 ) angle_rad_2d
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) p(dim_num)
  real ( kind = 8 ) p1(dim_num)
  real ( kind = 8 ) p2(dim_num)
  real ( kind = 8 ) p3(dim_num)

  p(1) = ( p1(1) - p2(1) ) * ( p3(1) - p2(1) ) &
       + ( p1(2) - p2(2) ) * ( p3(2) - p2(2) )

  p(2) =   ( p3(1) - p2(1) ) * ( p1(2) - p2(2) ) &
         - ( p1(1) - p2(1) ) * ( p3(2) - p2(2) ) 

  if ( p(1) == 0.0D+00 .and. p(2) == 0.0D+00 ) then

    angle_rad_2d = 0.0D+00

  else

    angle_rad_2d = atan2 ( p(2), p(1) )

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

  end if

  return
end
function atan4 ( y, x )

!*****************************************************************************80
!
!! ATAN4 computes the inverse tangent of the ratio Y / X.
!
!  Discussion:
!
!    ATAN4 returns an angle whose tangent is ( Y / X ), a job which
!    the built in functions ATAN and ATAN2 already do.
!
!    However:
!
!    * ATAN4 always returns a positive angle, between 0 and 2 PI,
!      while ATAN and ATAN2 return angles in the interval [-PI/2,+PI/2]
!      and [-PI,+PI] respectively;
!
!    * ATAN4 accounts for the signs of X and Y, (as does ATAN2).  The ATAN
!     function by contrast always returns an angle in the first or fourth
!     quadrants.
!
!  Modified:
!
!    15 October 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) Y, X, two quantities which represent the 
!    tangent of an angle.  If Y is not zero, then the tangent is (Y/X).
!
!    Output, real ( kind = 8 ) ATAN4, an angle between 0 and 2 * PI, whose
!    tangent is (Y/X), and which lies in the appropriate quadrant so that 
!    the signs of its cosine and sine match those of X and Y.
!
  implicit none

  real ( kind = 8 ) abs_x
  real ( kind = 8 ) abs_y
  real ( kind = 8 ) atan4
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) theta
  real ( kind = 8 ) theta_0
  real ( kind = 8 ) x
  real ( kind = 8 ) y
!
!  Special cases:
!
  if ( x == 0.0D+00 ) then

    if ( 0.0D+00 < y ) then
      theta = pi / 2.0D+00
    else if ( y < 0.0D+00 ) then
      theta = 3.0D+00 * pi / 2.0D+00
    else if ( y == 0.0D+00 ) then
      theta = 0.0D+00
    end if

  else if ( y == 0.0D+00 ) then

    if ( 0.0D+00 < x ) then
      theta = 0.0D+00
    else if ( x < 0.0D+00 ) then
      theta = pi
    end if
!
!  We assume that ATAN2 is correct when both arguments are positive.
!
  else

    abs_y = abs ( y )
    abs_x = abs ( x )

    theta_0 = atan2 ( abs_y, abs_x )

    if ( 0.0D+00 < x .and. 0.0D+00 < y ) then
      theta = theta_0
    else if ( x < 0.0D+00 .and. 0.0D+00 < y ) then
      theta = pi - theta_0
    else if ( x < 0.0D+00 .and. y < 0.0D+00 ) then
      theta = pi + theta_0
    else if ( 0.0D+00 < x .and. y < 0.0D+00 ) then
      theta = 2.0D+00 * pi - theta_0
    end if

  end if

  atan4 = theta

  return
end
function box_contains_point_2d ( box, p )

!*****************************************************************************80
!
!! BOX_CONTAINS_POINT_2D determines if a point is inside a box in 2D.
!
!  Discussion:
!
!    A unit box in 2D is the set of points (X,Y) with the property that
!
!      0.0 <= X <= 1.0
!    and
!      0.0 <= Y <= 1.0
!
!  Modified:
!
!    14 May 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) BOX(DIM_NUM,2), the lower left and upper right
!    corners of the box.
!
!    Input, real ( kind = 8 ) P(DIM_NUM), the point to be checked.
!
!    Output, logical BOX_CONTAINS_POINT_2D, is TRUE if the point is
!    inside the box.
!
  implicit none

  integer, parameter :: dim_num = 2

  real ( kind = 8 ) box(dim_num,2)
  logical box_contains_point_2d
  real ( kind = 8 ) p(dim_num)

  box_contains_point_2d = &
    all ( box(1:dim_num,1) <= p(1:dim_num) ) .and. &
    all (                     p(1:dim_num) <= box(1:dim_num,2) )

  return
end
subroutine circle_arc_point_near_2d ( r, c, theta1, theta2, x, nearest, &
  dist )

!*****************************************************************************80
!
!! CIRCLE_ARC_POINT_NEAR_2D : nearest point on a circular arc.
!
!  Discussion:
!
!    A circular arc is defined by the portion of a circle (R,C)
!    between two angles (THETA1,THETA2).
!
!    Thus, a point (X,Y) on a circular arc satisfies
!
!      ( X - C(1) ) * ( X - C(1) ) + ( Y - C(2) ) * ( Y - C(2) ) = R * R
!
!    and
!
!      Theta1 <= Theta <= Theta2
!
!  Modified:
!
!    22 October 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the circle.
!
!    Input, real ( kind = 8 ) C(2), the center of the circle.
!
!    Input, real ( kind = 8 ) THETA1, THETA2, the angles defining the arc,
!    in radians.  Normally, THETA1 < THETA2.
!
!    Input, real ( kind = 8 ) X(2), the point to be checked.
!
!    Output, real ( kind = 8 ) NEAREST(2), a point on the circular arc which is
!    nearest to the point.
!
!    Output, real ( kind = 8 ) DIST, the distance to the nearest point.
!
  implicit none

  integer, parameter :: dim_num = 2

  real ( kind = 8 ) atan4
  real ( kind = 8 ) c(dim_num)
  real ( kind = 8 ) dist
  real ( kind = 8 ) nearest(dim_num)
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) r
  real ( kind = 8 ) r2
  real ( kind = 8 ) r8_modp
  real ( kind = 8 ) theta
  real ( kind = 8 ) theta1
  real ( kind = 8 ) theta2
  real ( kind = 8 ) x(dim_num)
!
!  Special case, the zero circle.
!
  if ( r == 0.0D+00 ) then
    nearest(1:2) = c(1:2)
  end if
!
!  Determine the angle made by the point.
!
  theta = atan4 ( x(2) - c(2), x(1) - c(1) )
!
!  If the angle is between THETA1 and THETA2, then you can
!  simply project the point onto the arc.
!
  if ( r8_modp ( theta  - theta1,  2.0D+00 * pi ) <= &
       r8_modp ( theta2 - theta1,  2.0D+00 * pi ) ) then

    r2 = sqrt ( sum ( ( x(1:2) - c(1:2) )**2 ) )

    nearest(1:2) = c(1:2) + ( x(1:2) - c(1:2) ) * r / r2
!
!  Otherwise, if the angle is less than the negative of the
!  average of THETA1 and THETA2, it's on the side of the arc
!  where the endpoint associated with THETA2 is closest.
!
  else if ( r8_modp ( theta - 0.5D+00 * ( theta1 + theta2 ), 2.0D+00 * pi ) &
    <= pi ) then

    nearest(1:2) = c(1:2) + r * (/ cos ( theta2 ), sin ( theta2 ) /)
!
!  Otherwise, the endpoint associated with THETA1 is closest.
  else

    nearest(1:2) = c(1:2) + r * (/ cos ( theta1 ), sin ( theta1 ) /)

  end if

  dist = sqrt ( sum ( ( x(1:2) - nearest(1:2) )**2 ) )

  return
end
function circle_imp_contains_point_2d ( r, c, x )

!*****************************************************************************80
!
!! CIRCLE_IMP_CONTAINS_POINT_2D determines if a circle contains a point in 2D.
!
!  Discussion:
!
!    An implicit circle in 2D satisfies:
!
!      ( X - C(1) ) * ( X - C(1) ) + ( Y - C(2) ) * ( Y - C(2) ) = R * R
!
!  Modified:
!
!    19 October 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the circle.
!
!    Input, real ( kind = 8 ) C(2), the coordinates of the center of the circle.
!
!    Input, real ( kind = 8 ) X(2), the point to be checked.
!
!    Output, logical CIRCLE_IMP_CONTAINS_POINT_2D, is TRUE if the point
!    is inside or on the circle, FALSE otherwise.
!
  implicit none

  integer, parameter :: dim_num = 2

  real ( kind = 8 ) c(dim_num)
  logical circle_imp_contains_point_2d
  real ( kind = 8 ) r
  real ( kind = 8 ) x(dim_num)

  if ( ( x(1) - c(1) ) * ( x(1) - c(1) ) &
     + ( x(2) - c(2) ) * ( x(2) - c(2) ) <= r * r ) then
    circle_imp_contains_point_2d = .true.
  else
    circle_imp_contains_point_2d = .false.
  end if

  return
end
subroutine circle_imp_point_near_2d ( r, c, x, nearest, dist )

!*****************************************************************************80
!
!! CIRCLE_IMP_POINT_NEAR_2D: nearest ( implicit circle, point ) in 2D.
!
!  Discussion:
!
!    This routine finds the distance from a point to an implicitly
!    defined circle, and returns the point on the circle that is
!    nearest to the given point.
!
!    If the given point is the center of the circle, than any point
!    on the circle is "the" nearest.
!
!    An implicit circle in 2D satisfies the equation:
!
!      ( X - C(1) ) * ( X - C(1) ) + ( Y - C(2) ) * ( Y - C(2) ) = R * R
!
!  Modified:
!
!    19 October 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the circle.
!
!    Input, real ( kind = 8 ) C(2), the coordinates of the center of the circle.
!
!    Input, real ( kind = 8 ) X(2), the point to be checked.
!
!    Output, real ( kind = 8 ) NEAREST(2), the nearest point on the circle.
!
!    Output, real ( kind = 8 ) DIST, the distance of the point to the circle.
!
  implicit none

  integer, parameter :: dim_num = 2

  real ( kind = 8 ) c(dim_num)
  real ( kind = 8 ) dist
  real ( kind = 8 ) nearest(dim_num)
  real ( kind = 8 ) r
  real ( kind = 8 ) r2
  real ( kind = 8 ) x(dim_num)

  if ( all ( x(1:2) == c(1:2) ) ) then
    dist = r
    nearest(1:2) = c(1:2)
    nearest(1) = nearest(1) + r
    return
  end if

  r2 = sqrt ( sum ( ( x(1:2) - c(1:2) )**2 ) )

  dist = abs (  r2 - r )

  nearest(1:2) = c(1:2) + r * ( x(1:2) - c(1:2) ) / r2

  return
end
function circle_sector_contains_point_2d ( r, c, theta1, theta2, x )

!*****************************************************************************80
!
!! CIRCLE_SECTOR_CONTAINS_POINT_2D : is a point inside a circular sector?
!
!  Discussion:
!
!    A circular sector is formed by a circular arc, and the two straight line
!    segments that join its ends to the center of the circle.
!
!    A circular sector is defined by
!
!      ( X - C(1) ) * ( X - C(1) ) + ( Y - C(2) ) * ( Y - C(2) ) <= R * R
!
!    and
!
!      Theta1 <= Theta <= Theta2
!
!  Modified:
!
!    19 October 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the circle.
!
!    Input, real ( kind = 8 ) C(2), the coordinates of the center of the circle.
!
!    Input, real ( kind = 8 ) THETA1, THETA2, the angles defining the arc,
!    in radians.  Normally, THETA1 < THETA2.
!
!    Input, real ( kind = 8 ) X(2), the point to be checked.
!
!    Output, logical CIRCLE_SECTOR_CONTAINS_POINT_2D, is TRUE if the point
!    is inside or on the circular sector.
!
  implicit none

  integer, parameter :: dim_num = 2

  real ( kind = 8 ) atan4
  real ( kind = 8 ) c(dim_num)
  logical circle_sector_contains_point_2d
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) r
  real ( kind = 8 ) r8_modp
  real ( kind = 8 ) theta
  real ( kind = 8 ) theta1
  real ( kind = 8 ) theta2
  logical value
  real ( kind = 8 ) x(dim_num)

  value = .false.
!
!  Is the point inside the (full) circle?
!
  if ( ( x(1) - c(1) ) * ( x(1) - c(1) ) &
     + ( x(2) - c(2) ) * ( x(2) - c(2) ) <= r * r ) then
!
!  Is the point's angle within the arc's range?
!  Try to force the angles to lie between 0 and 2 * PI.
!
    theta = atan4 ( x(2) - c(2), x(1) - c(1) )

    if ( r8_modp ( theta  - theta1,  2.0D+00 * pi ) <= &
         r8_modp ( theta2 - theta1,  2.0D+00 * pi ) ) then

      value = .true.

    end if

  end if

  circle_sector_contains_point_2d = value

  return
end
function cos_deg ( angle )

!*****************************************************************************80
!
!! COS_DEG returns the cosine of an angle given in degrees.
!
!  Modified:
!
!    10 May 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ANGLE, the angle, in degrees.
!
!    Output, real ( kind = 8 ) COS_DEG, the cosine of the angle.
!
  implicit none

  real ( kind = 8 ) angle
  real ( kind = 8 ) cos_deg
  real ( kind = 8 ), parameter :: degrees_to_radians &
    = 3.141592653589793D+00 / 180.0D+00

  cos_deg = cos ( degrees_to_radians * angle )

  return
end
subroutine dtable_data_write ( output_unit, m, n, table )

!*****************************************************************************80
!
!! DTABLE_DATA_WRITE writes data to a double precision table file.
!
!  Discussion:
!
!    This routine writes a single line of output for each point,
!    containing its spatial coordinates.
!
!  Modified:
!
!    07 September 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer OUTPUT_UNIT, the output unit.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) TABLE(M,N), the table data.
!
  implicit none

  integer m
  integer n

  integer output_unit
  integer j
  character ( len = 30 ) string
  real ( kind = 8 ) table(m,n)
!
!  Create the format string.
!
  write ( string, '(a1,i8,a1,i8,a1,i8,a1)' ) '(', m, 'g', 14, '.', 6, ')'
  call s_blank_delete ( string )

  do j = 1, n
    write ( output_unit, string ) table(1:m,j)
  end do

  return
end
subroutine dtable_header_write ( output_file_name, output_unit, m, n )

!*****************************************************************************80
!
!! DTABLE_HEADER_WRITE writes the header to a double precision table file.
!
!  Modified:
!
!    07 September 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) OUTPUT_FILE_NAME, the output file name.
!
!    Input, integer OUTPUT_UNIT, the output unit.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
  implicit none

  integer m
  integer n
  character ( len = * ) output_file_name
  integer output_unit
  character ( len = 40 ) string
  real ( kind = 8 ), parameter :: x = 1.0D+00

  call timestring ( string )

  write ( output_unit, '(a)'       ) '#  ' // trim ( output_file_name )
  write ( output_unit, '(a)'       ) '#  created by TABLE_IO.F90'
  write ( output_unit, '(a)'       ) '#  at ' // trim ( string )
  write ( output_unit, '(a)'       ) '#'
  write ( output_unit, '(a,i8)'    ) '#  Spatial dimension M = ', m
  write ( output_unit, '(a,i8)'    ) '#  Number of points N = ', n
  write ( output_unit, '(a,g14.6)' ) '#  EPSILON (unit roundoff) = ', &
    epsilon ( x )
  write ( output_unit, '(a)'       ) '#'

  return
end
subroutine dtable_write ( output_file_name, m, n, table, header )

!*****************************************************************************80
!
!! DTABLE_WRITE writes a double precision table file.
!
!  Modified:
!
!    20 July 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) OUTPUT_FILE_NAME, the output file name.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) TABLE(M,N), the table data.
!
!    Input, logical HEADER, is TRUE if the header is to be included.
!
  implicit none

  integer m
  integer n

  logical header
  character ( len = * ) output_file_name
  integer output_status
  integer output_unit
  real ( kind = 8 ) table(m,n)

  call get_unit ( output_unit )

  open ( unit = output_unit, file = output_file_name, &
    status = 'replace', iostat = output_status )

  if ( output_status /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DTABLE_WRITE - Fatal error!'
    write ( *, '(a,i8)' ) '  Could not open the output file "' // &
      trim ( output_file_name ) // '" on unit ', output_unit
    stop
  end if

  if ( header ) then
    call dtable_header_write ( output_file_name, output_unit, m, n )
  end if

  call dtable_data_write ( output_unit, m, n, table )

  close ( unit = output_unit )

  return
end
subroutine file_name_inc ( file_name )

!*****************************************************************************80
!
!! FILE_NAME_INC generates the next filename in a series.
!
!  Discussion:
!
!    It is assumed that the digits in the name, whether scattered or
!    connected, represent a number that is to be increased by 1 on
!    each call.  If this number is all 9's on input, the output number
!    is all 0's.  Non-numeric letters of the name are unaffected, and
!    if the name contains no digits, then nothing is done.
!
!  Examples:
!
!      Input          Output
!      -----          ------
!      a7to11.txt     a7to12.txt
!      a7to99.txt     a8to00.txt
!      a9to99.txt     a0to00.txt
!      cat.txt        cat.txt
!
!  Modified:
!
!    12 October 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character ( len = * ) FILE_NAME.
!    On input, a character string to be incremented.
!    On output, the incremented string.
!
  implicit none

  character c
  integer digit
  character ( len = * ) file_name
  integer i
  integer lens

  lens = len_trim ( file_name )

  do i = lens, 1, -1

    c = file_name(i:i)

    if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then

      digit = ichar ( c ) - 48
      digit = digit + 1

      if ( digit == 10 ) then
        digit = 0
      end if

      c = char ( digit + 48 )

      file_name(i:i) = c

      if ( c /= '0' ) then
        return
      end if

    end if

  end do

  return
end
subroutine fmin_rc ( a, b, arg, status, value )

!*****************************************************************************80
!
!! FMIN_RC seeks a minimizer of a scalar function of a scalar variable.
!
!  Discussion:
!
!    FMIN_RC seeks an approximation to the point where F attains a minimum on
!    the interval (A,B).
!
!    The method used is a combination of golden section search and
!    successive parabolic interpolation.  Convergence is never much
!    slower than that for a Fibonacci search.  If F has a continuous
!    second derivative which is positive at the minimum (which is not
!    at A or B), then convergence is superlinear, and usually of the
!    order of about 1.324....
!
!    The routine is a revised version of the Brent FMIN algorithm,
!    which now uses reverse communication.
!
!    It is worth stating explicitly that this routine will NOT be
!    able to detect a minimizer that occurs at either initial endpoint
!    A or B.  If this is a concern to the user, then the user must
!    either ensure that the initial interval is larger, or to check
!    the function value at the returned minimizer against the values
!    at either endpoint.
!
!  Modified:
!
!    22 October 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization without Derivatives,
!    Prentice Hall, 1973.
!
!    David Kahaner, Clever Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1988.
!
!  Parameters
!
!    Input/output, real ( kind = 8 ) A, B.  On input, the left and right
!    endpoints of the initial interval.  On output, the lower and upper 
!    bounds for an interval containing the minimizer.  It is required 
!    that A < B.
!
!    Output, real ( kind = 8 ) ARG, the currently considered point.  The user
!    does not need to initialize this value.  On return with STATUS positive,
!    the user is requested to evaluate the function at ARG, and return
!    the value in VALUE.  On return with STATUS zero, ARG is the routine's
!    estimate for the function minimizer.
!
!    Input/output, integer STATUS, used to communicate between the user
!    and the routine.  The user only sets STATUS to zero on the first call,
!    to indicate that this is a startup call.  The routine returns STATUS
!    positive to request that the function be evaluated at ARG, or returns
!    STATUS as 0, to indicate that the iteration is complete and that
!    ARG is the estimated minimizer.
!
!    Input, real ( kind = 8 ) VALUE, the function value at ARG, as 
!    requested by the routine on the previous call.
!
!  Local parameters:
!
!    C is the squared inverse of the golden ratio.
!
!    EPS is the square root of the relative machine precision.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) arg
  real ( kind = 8 ) b
  real ( kind = 8 ), save :: c
  real ( kind = 8 ), save :: d
  real ( kind = 8 ), save :: e
  real ( kind = 8 ), save :: eps
  real ( kind = 8 ), save :: fu
  real ( kind = 8 ), save :: fv
  real ( kind = 8 ), save :: fw
  real ( kind = 8 ), save :: fx
  real ( kind = 8 ), save :: midpoint
  real ( kind = 8 ), save :: p
  real ( kind = 8 ), save :: q
  real ( kind = 8 ), save :: r
  integer status
  real ( kind = 8 ), save :: tol
  real ( kind = 8 ), save :: tol1
  real ( kind = 8 ), save :: tol2
  real ( kind = 8 ), save :: u
  real ( kind = 8 ), save :: v
  real ( kind = 8 ) value
  real ( kind = 8 ), save :: w
  real ( kind = 8 ), save :: x
!
!  STATUS (INPUT) = 0, startup.
!
  if ( status == 0 ) then

    if ( b <= a ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'FMIN_RC - Fatal error!'
      write ( *, '(a)' ) '  A < B is required, but'
      write ( *, '(a,g14.6)' ) '  A = ', a
      write ( *, '(a,g14.6)' ) '  B = ', b
      status = -1
      stop
    end if

    c = 0.5D+00 * ( 3.0D+00 - sqrt ( 5.0D+00 ) )

    eps = sqrt ( epsilon ( eps ) )
    tol = epsilon ( tol )

    v = a + c * ( b - a )
    w = v
    x = v
    e = 0.0D+00

    status = 1
    arg = x

    return
!
!  STATUS (INPUT) = 1, return with initial function value of FX.
!
  else if ( status == 1 ) then

    fx = value
    fv = fx
    fw = fx
!
!  STATUS (INPUT) = 2 or more, update the data.
!
  else if ( 2 <= status ) then

    fu = value

    if ( fu <= fx ) then

      if ( x <= u ) then
        a = x
      else
        b = x
      end if

      v = w
      fv = fw
      w = x
      fw = fx
      x = u
      fx = fu

    else

      if ( u < x ) then
        a = u
      else
        b = u
      end if

      if ( fu <= fw .or. w == x ) then
        v = w
        fv = fw
        w = u
        fw = fu
      else if ( fu <= fv .or. v == x .or. v == w ) then
        v = u
        fv = fu
      end if

    end if

  end if
!
!  Take the next step.
!
  midpoint = 0.5D+00 * ( a + b )
  tol1 = eps * abs ( x ) + tol / 3.0D+00
  tol2 = 2.0D+00 * tol1
!
!  If the stopping criterion is satisfied, we can exit.
!
  if ( abs ( x - midpoint ) <= ( tol2 - 0.5D+00 * ( b - a ) ) ) then
    status = 0
    return
  end if
!
!  Is golden-section necessary?
!
  if ( abs ( e ) <= tol1 ) then
    if ( midpoint <= x ) then
      e = a - x
    else
      e = b - x
    end if

    d = c * e
!
!  Consider fitting a parabola.
!
  else

    r = ( x - w ) * ( fx - fv )
    q = ( x - v ) * ( fx - fw )
    p = ( x - v ) * q - ( x - w ) * r
    q = 2.0D+00 * ( q - r )
    if ( 0.0D+00 < q ) then
      p = -p
    end if
    q = abs ( q )
    r = e
    e = d
!
!  Choose a golden-section step if the parabola is not advised.
!
    if ( &
      ( abs ( 0.5D+00 * q * r ) <= abs ( p ) ) .or. &
      ( p <= q * ( a - x ) ) .or. &
      ( q * ( b - x ) <= p ) ) then

      if ( midpoint <= x ) then
        e = a - x
      else
        e = b - x
      end if

      d = c * e
!
!  Choose a parabolic interpolation step.
!
    else

      d = p / q
      u = x + d

      if ( ( u - a ) < tol2 ) then
        d = sign ( tol1, midpoint - x )
      end if

      if ( ( b - u ) < tol2 ) then
        d = sign ( tol1, midpoint - x )
      end if

    end if

  end if
!
!  F must not be evaluated too close to X.
!
  if ( tol1 <= abs ( d ) ) then
    u = x + d
  end if

  if ( abs ( d ) < tol1 ) then
    u = x + sign ( tol1, d )
  end if
!
!  Request value of F(U).
!
  arg = u
  status = status + 1

  return
end
subroutine get_unit ( iunit )

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

  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

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

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

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

    end if

  end do

  return
end
subroutine hex_grid_angle ( box, center, angle, h, n, r )

!*****************************************************************************80
!
!! HEX_GRID_ANGLE sets the points in an angled hex grid in a box.
!
!  Discussion:
!
!    DIM_NUM = 2
!
!  Modified:
!
!    14 May 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) BOX(DIM_NUM,2), the lower left and upper right
!    corners of the box.
!
!    Input, real ( kind = 8 ) CENTER(DIM_NUM), the center of the grid.
!    This point must be inside the unit square.
!
!    Input, real ( kind = 8 ) ANGLE, the angle, in degrees, of the grid.
!    Normally, 0 <= ANGLE <= 180, but any value is allowed.
!
!    Input, real ( kind = 8 ) H, the spacing between neighboring
!    points on a grid line.
!
!    Input, integer N, the number of points of the angled hex grid
!    that are within the unit square.  This value may have been computed
!    by calling HEX_GRID_ANGLE_01_SIZE.
!
!    Output, real ( kind = 8 ) R(DIM_NUM,N), the grid points.
!
  implicit none

  integer n
  integer, parameter :: dim_num = 2

  real ( kind = 8 ) angle
  real ( kind = 8 ) angle2
  real ( kind = 8 ) box(dim_num,2)
  logical box_contains_point_2d
  real ( kind = 8 ) center(dim_num)
  real ( kind = 8 ) cos_deg
  real ( kind = 8 ) h
  integer i
  integer j
  integer k
  integer layer
  integer layer_size
  real ( kind = 8 ) point(dim_num)
  real ( kind = 8 ) r(dim_num,n)
  real ( kind = 8 ) r8_modp
  real ( kind = 8 ) sin_deg
  integer size
!
!  Ninny checks.
!
  if ( .not. box_contains_point_2d ( box, center ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HEX_GRID_ANGLE - Fatal error!'
    write ( *, '(a)' ) '  The center point of the grid is not'
    write ( *, '(a)' ) '  inside the box.'
    write ( *, '(a,2g14.6)' ) '  CENTER = ', center(1:dim_num)
    stop
  end if

  if ( h == 0.0D+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HEX_GRID_ANGLE - Fatal error!'
    write ( *, '(a)' ) '  The grid spacing must be nonzero.'
    write ( *, '(a,g14.6)' ) '  H = ', h
    stop
  end if

  layer = 0
  point(1:dim_num) = center(1:dim_num)

  k = 1
  if ( k <= n ) then
    r(1:dim_num,k) = center(1:dim_num)
  end if

  do

    layer = layer + 1

    layer_size = 0

    angle2 = angle
!
!  Compute the first point on the new layer.
!
    point(1:dim_num) = point(1:dim_num) &
      + h * (/ cos_deg ( angle2 ), sin_deg ( angle2 ) /)

    angle2 = r8_modp ( angle2 + 60.0D+00, 360.0D+00 )

    do i = 1, 6

      angle2 = r8_modp ( angle2 + 60.0D+00, 360.0D+00 )

      do j = 1, layer

        point(1:dim_num) = point(1:dim_num) &
          + h * (/ cos_deg ( angle2 ), sin_deg ( angle2 ) /)

        if ( box_contains_point_2d ( box, point ) ) then

          layer_size = layer_size + 1
          k = k + 1

          if ( k <= n ) then
            r(1:dim_num,k) = point(1:dim_num)
          end if

        end if

      end do

    end do

    if ( layer_size == 0 ) then
      exit
    end if

  end do

  return
end
subroutine hex_grid_angle_size ( box, center, angle, h, n )

!*****************************************************************************80
!
!! HEX_GRID_ANGLE_SIZE counts the points in an angled hex grid in a box.
!
!  Discussion:
!
!    DIM_NUM = 2
!
!  Modified:
!
!    14 May 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) BOX(DIM_NUM,2), the lower left and upper right
!    corners of the box.
!
!    Input, real ( kind = 8 ) CENTER(DIM_NUM), the center of the grid.
!    This point must be inside the box
!
!    Input, real ( kind = 8 ) ANGLE, the angle, in degrees, of the grid.
!    Normally, 0 <= ANGLE <= 180, but any value is allowed.
!
!    Input, real ( kind = 8 ) H, the spacing between neighboring
!    points on a grid line.
!
!    Output, integer N, the number of points of the angled hex grid
!    that are within the unit square.
!
  implicit none

  integer, parameter :: dim_num = 2

  real ( kind = 8 ) angle
  real ( kind = 8 ) angle2
  real ( kind = 8 ) box(dim_num,2)
  logical box_contains_point_2d
  real ( kind = 8 ) center(dim_num)
  real ( kind = 8 ) cos_deg
  real ( kind = 8 ) h
  integer i
  integer j
  integer layer
  integer layer_size
  integer n
  real ( kind = 8 ) point(dim_num)
  real ( kind = 8 ) r8_modp
  real ( kind = 8 ) sin_deg
!
!  Ninny checks.
!
  if ( .not. box_contains_point_2d ( box, center ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HEX_GRID_ANGLE_SIZE - Fatal error!'
    write ( *, '(a)' ) '  The center point of the grid is not'
    write ( *, '(a)' ) '  inside the box.'
    write ( *, '(a,2g14.6)' ) '  CENTER = ', center(1:dim_num)
    stop
  end if

  if ( h == 0.0D+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HEX_GRID_ANGLE_SIZE - Fatal error!'
    write ( *, '(a)' ) '  The grid spacing must be nonzero.'
    write ( *, '(a,g14.6)' ) '  H = ', h
    stop
  end if

  n = 0

  layer = 0
  point(1:dim_num) = center(1:dim_num)

  n = 1

  do

    layer = layer + 1

    layer_size = 0

    angle2 = angle
!
!  Compute the first point on the new layer.
!
    point(1:dim_num) = point(1:dim_num) &
      + h * (/ cos_deg ( angle2 ), sin_deg ( angle2 ) /)

    angle2 = r8_modp ( angle2 + 60.0D+00, 360.0D+00 )

    do i = 1, 6

      angle2 = r8_modp ( angle2 + 60.0D+00, 360.0D+00 )

      do j = 1, layer

        point(1:dim_num) = point(1:dim_num) &
          + h * (/ cos_deg ( angle2 ), sin_deg ( angle2 ) /)

        if ( box_contains_point_2d ( box, point ) ) then
          layer_size = layer_size + 1
          n = n + 1
        end if

      end do

    end do

    if ( layer_size == 0 ) then
      exit
    end if

  end do

  return
end
function hexagon_contains_point_2d ( v, p )

!*****************************************************************************80
!
!! HEXAGON_CONTAINS_POINT_2D finds if a point is inside a hexagon in 2D.
!
!  Discussion:
!
!    This test is only valid if the hexagon is convex.
!
!  Modified:
!
!    20 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) V(2,6), the vertics, in counterclockwise order.
!
!    Input, real ( kind = 8 ) P(2), the point to be tested.
!
!    Output, logical HEXAGON_CONTAINS_POINT_2D, is TRUE if X is in the hexagon.
!
  implicit none

  integer, parameter :: dim_num = 2
  integer, parameter :: n = 6

  logical hexagon_contains_point_2d
  integer i
  integer j
  real ( kind = 8 ) p(dim_num)
  real ( kind = 8 ) v(dim_num,n)
!
!  A point is inside a convex hexagon if and only if it is "inside"
!  each of the 6 halfplanes defined by lines through consecutive
!  vertices.
!
  do i = 1, n

    j = mod ( i, n ) + 1

    if (  v(1,i) * ( v(2,j) - p(2  ) ) &
        + v(1,j) * ( p(2  ) - v(2,i) ) &
        + p(1  ) * ( v(2,i) - v(2,j) ) < 0.0D+00 ) then

      hexagon_contains_point_2d = .false.
      return

    end if

  end do

  hexagon_contains_point_2d = .true.

  return
end
function i4_modp ( i, j )

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

  integer i
  integer i4_modp
  integer j

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

  i4_modp = mod ( i, j )

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

  return
end
function i4_wrap ( ival, ilo, ihi )

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

  integer i4_modp
  integer i4_wrap
  integer ihi
  integer ilo
  integer ival
  integer jhi
  integer jlo
  integer wide

  jlo = min ( ilo, ihi )
  jhi = max ( ilo, ihi )

  wide = jhi - jlo + 1

  if ( wide == 1 ) then
    i4_wrap = jlo
  else
    i4_wrap = jlo + i4_modp ( ival - jlo, wide )
  end if

  return
end
subroutine p00_boundary_eps ( test, h, show_nodes, eps_file_name )

!*****************************************************************************80
!
!! P00_BOUNDARY_EPS draws the boundary of a region as an EPS file.
!
!  Modified:
!
!    20 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the number of the problem.
!
!    Input, real ( kind = 8 ) H, the maximum size of a segment of the boundary.
!    This controls how smoothly curved sections of the boundary will be drawn.
!
!    Input, logical SHOW_NODES, is TRUE if the boundary nodes are
!    to be shown.
!
!    Input, character ( len = * ) EPS_FILE_NAME, the name of the EPS file
!    to create.
!
  implicit none

  integer, parameter :: m = 2

  real ( kind = 8 ), allocatable, dimension ( :, : ) :: boundary
  character ( len = * ) eps_file_name
  integer eps_unit
  real ( kind = 8 ) h
  real ( kind = 8 ), allocatable, dimension ( : ) :: hi
  integer ios
  integer j
  real ( kind = 8 ), allocatable, dimension ( : ) :: lo
  integer pass
  real ( kind = 8 ) scale
  integer segment
  integer segment_length
  integer segment_num
  logical show_nodes
  character ( len = 40 ) string
  integer test
  real ( kind = 8 ) x_max
  real ( kind = 8 ) x_min
  integer x_ps
  integer :: x_ps_max = 576
  integer :: x_ps_max_clip = 594
  integer :: x_ps_max_user  
  integer :: x_ps_min = 36
  integer :: x_ps_min_clip = 18
  integer :: x_ps_min_user  
  real ( kind = 8 ) x_scale
  real ( kind = 8 ) y_max
  real ( kind = 8 ) y_min
  integer y_ps
  integer :: y_ps_max = 666
  integer :: y_ps_max_clip = 684
  integer :: y_ps_max_user  
  integer :: y_ps_min = 126
  integer :: y_ps_min_clip = 108
  integer :: y_ps_min_user
  real ( kind = 8 ) y_scale

  call p00_boundary_segment_num ( test, segment_num )

  allocate ( lo(1:m) )
  allocate ( hi(1:m) )

  call p00_box ( test, m, lo, hi )

  x_min = lo(1) - 0.025D+00 * ( hi(1) - lo(1) )
  y_min = lo(2) - 0.025D+00 * ( hi(2) - lo(2) )
  x_max = hi(1) + 0.025D+00 * ( hi(1) - lo(1) )
  y_max = hi(2) + 0.025D+00 * ( hi(2) - lo(2) )

  x_scale = x_max - x_min
  y_scale = y_max - y_min
  scale = max ( x_scale, y_scale )
!
!  Determine the PostScript coordinates of the used box.
!
  x_ps_min_user = ( ( x_ps_max + x_ps_min ) &
    - int ( real ( x_ps_max - x_ps_min, kind = 8 ) * x_scale / scale ) ) / 2
  x_ps_max_user = ( ( x_ps_max + x_ps_min ) &
    + int ( real ( x_ps_max - x_ps_min, kind = 8 ) * x_scale / scale ) ) / 2
  y_ps_min_user = ( ( y_ps_max + y_ps_min ) &
    - int ( real ( y_ps_max - y_ps_min, kind = 8 ) * y_scale / scale ) ) / 2
  y_ps_max_user = ( ( y_ps_max + y_ps_min ) &
    + int ( real ( y_ps_max - y_ps_min, kind = 8 ) * y_scale / scale ) ) / 2

  if ( x_scale < y_scale ) then
    x_max = x_max + 0.5D+00 * ( y_scale - x_scale )
    x_min = x_min - 0.5D+00 * ( y_scale - x_scale )
  else if ( y_scale < x_scale ) then
    y_max = y_max + 0.5D+00 * ( x_scale - y_scale )
    y_min = y_min - 0.5D+00 * ( x_scale - y_scale )
  end if

  call get_unit ( eps_unit )

  open ( unit = eps_unit, file = eps_file_name, status = 'replace', &
    iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_BOUNDARY_EPS - Fatal error!'
    write ( *, '(a)' ) '  Could not open the output EPS file:'
    write ( *, '(a)' ) '    "' // trim ( eps_file_name ) // '".'
    write ( *, '(a,i12)' ) '  IOS = ', ios
    stop
  end if

  call timestring ( string )

  write ( eps_unit, '(a)' ) '%!PS-Adobe-3.0 EPSF-3.0'
  write ( eps_unit, '(a)' ) '%%Creator: p00_boundary_eps.f90'
  write ( eps_unit, '(a)' ) '%%Title: ' // trim ( eps_file_name )
  write ( eps_unit, '(a)' ) '%%CreationDate: ' // trim ( string )
  write ( eps_unit, '(a)' ) '%%Pages: 1'
  write ( eps_unit, '(a,i4,2x,i4,2x,i4,2x,i4)' ) '%%BoundingBox: ', &
    x_ps_min_user, y_ps_min_user, x_ps_max_user, y_ps_max_user
  write ( eps_unit, '(a)' ) '%%Document-Fonts: Times-Roman'
  write ( eps_unit, '(a)' ) '%%LanguageLevel: 1'
  write ( eps_unit, '(a)' ) '%%EndComments'
  write ( eps_unit, '(a)' ) '%%BeginProlog'
  write ( eps_unit, '(a)' ) '/inch {72 mul} def'
  write ( eps_unit, '(a)' ) '%%EndProlog'
  write ( eps_unit, '(a)' ) '%%Page: 1 1'
  write ( eps_unit, '(a)' ) 'save'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Set the RGB line color to very light gray.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '0.900  0.900  0.900 setrgbcolor'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Draw a gray border around the user box.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) 'newpath'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_min_user, y_ps_min_user, ' moveto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_max_user, y_ps_min_user, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_max_user, y_ps_max_user, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_min_user, y_ps_max_user, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_min_user, y_ps_min_user, ' lineto'
  write ( eps_unit, '(a)' ) 'stroke'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Set the RGB line color to black.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '0.000  0.000  0.000 setrgbcolor'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Set the font and its size.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '/Times-Roman findfont'
  write ( eps_unit, '(a)' ) '0.50 inch scalefont'
  write ( eps_unit, '(a)' ) 'setfont'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Define a clipping polygon.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) 'newpath'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_min_clip, y_ps_min_clip, ' moveto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_max_clip, y_ps_min_clip, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_max_clip, y_ps_max_clip, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_min_clip, y_ps_max_clip, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_min_clip, y_ps_min_clip, ' lineto'
  write ( eps_unit, '(a)' ) 'clip newpath'

  do segment = 1, segment_num

    call p00_boundary_segment_length ( test, segment, h, segment_length )

    allocate ( boundary(1:m,1:segment_length) )

    call p00_boundary_segment ( test, segment, m, segment_length, boundary )

    if ( show_nodes ) then

      write ( eps_unit, '(a)' ) '%'
      write ( eps_unit, '(a)' ) '%  Set the RGB line color to green.'
      write ( eps_unit, '(a)' ) '%'
      write ( eps_unit, '(a)' ) '0.000  0.750  0.150 setrgbcolor'
      write ( eps_unit, '(a)' ) '%'
      write ( eps_unit, '(a)' ) '%  Draw the nodes.'
      write ( eps_unit, '(a)' ) '%'

      do j = 1, segment_length

        x_ps = int ( &
          ( ( x_max - boundary(1,j)         ) * real ( x_ps_min, kind = 8 )   &
          + (         boundary(1,j) - x_min ) * real ( x_ps_max, kind = 8 ) ) &
          / ( x_max                 - x_min ) )

        y_ps = int ( &
          ( ( y_max - boundary(2,j)         ) * real ( y_ps_min, kind = 8 )   &
          + (         boundary(2,j) - y_min ) * real ( y_ps_max, kind = 8 ) ) &
          / ( y_max                 - y_min ) )

        write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) 'newpath ', x_ps, y_ps, &
          ' 5 0 360 arc closepath fill'

      end do
    end if

    write ( eps_unit, '(a)' ) '%'
    write ( eps_unit, '(a)' ) '%  Set the RGB line color to red.'
    write ( eps_unit, '(a)' ) '%'
    write ( eps_unit, '(a)' ) '0.900  0.200  0.100 setrgbcolor'
    write ( eps_unit, '(a)' ) '%'
    write ( eps_unit, '(a)' ) '%  Increase the linewidth to 3.'
    write ( eps_unit, '(a)' ) '%'
    write ( eps_unit, '(a)' ) '3  setlinewidth'
    write ( eps_unit, '(a)' ) '%'
    write ( eps_unit, '(a)' ) '%  Draw the boundary lines.'
    write ( eps_unit, '(a)' ) '%'

    do j = 1, segment_length
      x_ps = int ( &
        ( ( x_max - boundary(1,j)         ) * real ( x_ps_min, kind = 8 )   &
        + (         boundary(1,j) - x_min ) * real ( x_ps_max, kind = 8 ) ) &
        / ( x_max                 - x_min ) )

      y_ps = int ( &
        ( ( y_max - boundary(2,j)         ) * real ( y_ps_min, kind = 8 )   &
        + (         boundary(2,j) - y_min ) * real ( y_ps_max, kind = 8 ) ) &
        / ( y_max                 - y_min ) )

      if ( j == 1 ) then
        write ( eps_unit, '(i4,2x,i4,2x,a)' ) x_ps, y_ps, ' moveto'
      else
        write ( eps_unit, '(i4,2x,i4,2x,a)' ) x_ps, y_ps, ' lineto'
      end if

    end do

    write ( eps_unit, '(a)' ) 'stroke'

    deallocate ( boundary )

  end do

  deallocate ( hi )
  deallocate ( lo )

  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) 'restore showpage'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '% End of page'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%%Trailer'
  write ( eps_unit, '(a)' ) '%%EOF'

  close ( unit = eps_unit )

  if ( .false. ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_BOUNDARY_EPS:'
    write ( *, '(a)' ) '  An encapsulated PostScript file was created'
    write ( *, '(a)' ) '  containing an image of the boundary.'
    write ( *, '(a)' ) '  The file is named "' // trim ( eps_file_name ) // '".'
  end if

  return
end
subroutine p00_boundary_nearest ( test, dim_num, n, point, boundary )

!*****************************************************************************80
!
!! P00_BOUNDARY_NEAREST returns a nearest boundary point for any problem.
!
!  Discussion:
!
!    The given input point need not be inside the region.
!
!    In some cases, more than one boundary point may be "nearest",
!    but only one will be returned.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(DIM_NUM,N), the coordinates 
!    of the points.
!
!    Output, real ( kind = 8 ) BOUNDARY(DIM_NUM,N), points on the boundary
!    that are nearest to each point.
!
  implicit none

  integer dim_num
  integer n

  real ( kind = 8 ), dimension ( dim_num, n ) :: boundary
  real ( kind = 8 ), dimension ( dim_num, n ) :: point
  integer test

  if ( test == 1 ) then
    call p01_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 2 ) then
    call p02_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 3 ) then
    call p03_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 4 ) then
    call p04_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 5 ) then
    call p05_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 6 ) then
    call p06_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 7 ) then
    call p07_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 8 ) then
    call p08_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 9 ) then
    call p09_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 10 ) then
    call p10_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 11 ) then
    call p11_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 12 ) then
    call p12_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 13 ) then
    call p13_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 14 ) then
    call p14_boundary_nearest ( dim_num, n, point, boundary )
  else if ( test == 15 ) then
    call p15_boundary_nearest ( dim_num, n, point, boundary )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_BOUNDARY_NEAREST - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_boundary_project ( test, dim_num, n, point )

!*****************************************************************************80
!
!! P00_BOUNDARY_PROJECT projects exterior points to the boundary.
!
!  Discussion:
!
!    Interior points are not changed.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input/output, real ( kind = 8 ) POINT(DIM_NUM,N), the coordinates 
!    of the points.  Any input points that are exterior to the region
!    are replaced on output by the nearest boundary point.
!
  implicit none

  integer dim_num
  integer n

  real ( kind = 8 ), dimension ( dim_num, n ) :: point
  integer test

  if ( test == 1 ) then
    call p01_boundary_project ( dim_num, n, point )
  else if ( test == 2 ) then
    call p02_boundary_project ( dim_num, n, point )
  else if ( test == 3 ) then
    call p03_boundary_project ( dim_num, n, point )
  else if ( test == 4 ) then
    call p04_boundary_project ( dim_num, n, point )
  else if ( test == 5 ) then
    call p05_boundary_project ( dim_num, n, point )
  else if ( test == 6 ) then
    call p06_boundary_project ( dim_num, n, point )
  else if ( test == 7 ) then
    call p07_boundary_project ( dim_num, n, point )
  else if ( test == 8 ) then
    call p08_boundary_project ( dim_num, n, point )
  else if ( test == 9 ) then
    call p09_boundary_project ( dim_num, n, point )
  else if ( test == 10 ) then
    call p10_boundary_project ( dim_num, n, point )
  else if ( test == 11 ) then
    call p11_boundary_project ( dim_num, n, point )
  else if ( test == 12 ) then
    call p12_boundary_project ( dim_num, n, point )
  else if ( test == 13 ) then
    call p13_boundary_project ( dim_num, n, point )
  else if ( test == 14 ) then
    call p14_boundary_project ( dim_num, n, point )
  else if ( test == 15 ) then
    call p15_boundary_project ( dim_num, n, point )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_BOUNDARY_PROJECT - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_boundary_segment ( test, segment_index, m, &
  segment_length, segment )

!*****************************************************************************80
!
!! P00_BOUNDARY_SEGMENT returns a boundary segment in a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer SEGMENT_INDEX, the index of the segment.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer SEGMENT_LENGTH, the number of points in the segment.
!
!    Output, real ( kind = 8 ) SEGMENT(M,SEGMENT_LENGTH), the 
!    points that make up the boundary segment.
!
  implicit none

  integer m
  integer segment_length

  integer segment_index
  real ( kind = 8 ) segment(m,segment_length)
  integer test

  if ( test == 1 ) then
    call p01_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 2 ) then
    call p02_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 3 ) then
    call p03_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 4 ) then
    call p04_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 5 ) then
    call p05_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 6 ) then
    call p06_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 7 ) then
    call p07_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 8 ) then
    call p08_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 9 ) then
    call p09_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 10 ) then
    call p10_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 11 ) then
    call p11_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 12 ) then
    call p12_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 13 ) then
    call p13_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 14 ) then
    call p14_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else if ( test == 15 ) then
    call p15_boundary_segment ( segment_index, m, segment_length, &
      segment )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_BOUNDARY_SEGMENT - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_boundary_segment_length ( test, segment_index, h, &
  segment_length )

!*****************************************************************************80
!
!! P00_BOUNDARY_SEGMENT_LENGTH returns boundary segment lengths in a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer SEGMENT_INDEX, the index of one of the boundary segments.
!
!    Input, real ( kind = 8 ) H, the suggested spacing between points.
!
!    Output, integer SEGMENT_LENGTH, the number of points in the segment.
!
  implicit none

  real ( kind = 8 ) h
  integer segment_index
  integer segment_length
  integer test

  if ( test == 1 ) then
    call p01_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 2 ) then
    call p02_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 3 ) then
    call p03_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 4 ) then
    call p04_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 5 ) then
    call p05_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 6 ) then
    call p06_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 7 ) then
    call p07_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 8 ) then
    call p08_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 9 ) then
    call p09_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 10 ) then
    call p10_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 11 ) then
    call p11_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 12 ) then
    call p12_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 13 ) then
    call p13_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 14 ) then
    call p14_boundary_segment_length ( segment_index, h, segment_length )
  else if ( test == 15 ) then
    call p15_boundary_segment_length ( segment_index, h, segment_length )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_BOUNDARY_SEGMENT_LENGTH - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_boundary_segment_num ( test, boundary_segment_num )

!*****************************************************************************80
!
!! P00_BOUNDARY_SEGMENT_NUM counts the boundary segments in a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Output, integer BOUNDARY_SEGMENT_NUM, the number of boundary segments.
!
  implicit none

  integer boundary_segment_num
  integer test

  if ( test == 1 ) then
    call p01_boundary_segment_num ( boundary_segment_num )
  else if ( test == 2 ) then
    call p02_boundary_segment_num ( boundary_segment_num )
  else if ( test == 3 ) then
    call p03_boundary_segment_num ( boundary_segment_num )
  else if ( test == 4 ) then
    call p04_boundary_segment_num ( boundary_segment_num )
  else if ( test == 5 ) then
    call p05_boundary_segment_num ( boundary_segment_num )
  else if ( test == 6 ) then
    call p06_boundary_segment_num ( boundary_segment_num )
  else if ( test == 7 ) then
    call p07_boundary_segment_num ( boundary_segment_num )
  else if ( test == 8 ) then
    call p08_boundary_segment_num ( boundary_segment_num )
  else if ( test == 9 ) then
    call p09_boundary_segment_num ( boundary_segment_num )
  else if ( test == 10 ) then
    call p10_boundary_segment_num ( boundary_segment_num )
  else if ( test == 11 ) then
    call p11_boundary_segment_num ( boundary_segment_num )
  else if ( test == 12 ) then
    call p12_boundary_segment_num ( boundary_segment_num )
  else if ( test == 13 ) then
    call p13_boundary_segment_num ( boundary_segment_num )
  else if ( test == 14 ) then
    call p14_boundary_segment_num ( boundary_segment_num )
  else if ( test == 15 ) then
    call p15_boundary_segment_num ( boundary_segment_num )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_BOUNDARY_SEGMENT_NUM - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_box ( test, m, lo, hi )

!*****************************************************************************80
!
!! P00_BOX returns a bounding box for a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer M, the spatial dimension.
!
!    Output, real ( kind = 8 ) LO(M), HI(M), the lower and
!    upper corners of a bounding box.
!
  implicit none

  integer m

  real ( kind = 8 ) hi(m)
  real ( kind = 8 ) lo(m)
  integer test

  if ( test == 1 ) then
    call p01_box ( m, lo, hi )
  else if ( test == 2 ) then
    call p02_box ( m, lo, hi )
  else if ( test == 3 ) then
    call p03_box ( m, lo, hi )
  else if ( test == 4 ) then
    call p04_box ( m, lo, hi )
  else if ( test == 5 ) then
    call p05_box ( m, lo, hi )
  else if ( test == 6 ) then
    call p06_box ( m, lo, hi )
  else if ( test == 7 ) then
    call p07_box ( m, lo, hi )
  else if ( test == 8 ) then
    call p08_box ( m, lo, hi )
  else if ( test == 9 ) then
    call p09_box ( m, lo, hi )
  else if ( test == 10 ) then
    call p10_box ( m, lo, hi )
  else if ( test == 11 ) then
    call p11_box ( m, lo, hi )
  else if ( test == 12 ) then
    call p12_box ( m, lo, hi )
  else if ( test == 13 ) then
    call p13_box ( m, lo, hi )
  else if ( test == 14 ) then
    call p14_box ( m, lo, hi )
  else if ( test == 15 ) then
    call p15_box ( m, lo, hi )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_BOX - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_density ( test, m, n, point, density )

!*****************************************************************************80
!
!! P00_DENSITY returns the density for a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
!    Output, real ( kind = 8 ) DENSITY(N), the mesh density 
!    at each point.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ) density(n)
  real ( kind = 8 ) point(m,n)
  integer test

  if ( test == 1 ) then
    call p01_density ( m, n, point, density )
  else if ( test == 2 ) then
    call p02_density ( m, n, point, density )
  else if ( test == 3 ) then
    call p03_density ( m, n, point, density )
  else if ( test == 4 ) then
    call p04_density ( m, n, point, density )
  else if ( test == 5 ) then
    call p05_density ( m, n, point, density )
  else if ( test == 6 ) then
    call p06_density ( m, n, point, density )
  else if ( test == 7 ) then
    call p07_density ( m, n, point, density )
  else if ( test == 8 ) then
    call p08_density ( m, n, point, density )
  else if ( test == 9 ) then
    call p09_density ( m, n, point, density )
  else if ( test == 10 ) then
    call p10_density ( m, n, point, density )
  else if ( test == 11 ) then
    call p11_density ( m, n, point, density )
  else if ( test == 12 ) then
    call p12_density ( m, n, point, density )
  else if ( test == 13 ) then
    call p13_density ( m, n, point, density )
  else if ( test == 14 ) then
    call p14_density ( m, n, point, density )
  else if ( test == 15 ) then
    call p15_density ( m, n, point, density )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_DENSITY - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_element_size ( test, element_size )

!*****************************************************************************80
!
!! P00_ELEMENT_SIZE returns a typical element size for a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, real ( kind = 8 ) ELEMENT_SIZE, a typical element size.
!
  implicit none

  real ( kind = 8 ) element_size
  integer test

  if ( test == 1 ) then
    call p01_element_size ( element_size )
  else if ( test == 2 ) then
    call p02_element_size ( element_size )
  else if ( test == 3 ) then
    call p03_element_size ( element_size )
  else if ( test == 4 ) then
    call p04_element_size ( element_size )
  else if ( test == 5 ) then
    call p05_element_size ( element_size )
  else if ( test == 6 ) then
    call p06_element_size ( element_size )
  else if ( test == 7 ) then
    call p07_element_size ( element_size )
  else if ( test == 8 ) then
    call p08_element_size ( element_size )
  else if ( test == 9 ) then
    call p09_element_size ( element_size )
  else if ( test == 10 ) then
    call p10_element_size ( element_size )
  else if ( test == 11 ) then
    call p11_element_size ( element_size )
  else if ( test == 12 ) then
    call p12_element_size ( element_size )
  else if ( test == 13 ) then
    call p13_element_size ( element_size )
  else if ( test == 14 ) then
    call p14_element_size ( element_size )
  else if ( test == 15 ) then
    call p15_element_size ( element_size )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_ELEMENT_SIZE - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_fixed_num ( test, fixed_num )

!*****************************************************************************80
!
!! P00_FIXED_NUM returns the number of fixed points in a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Output, integer FIXED_NUM, the number of fixed points.
!
  implicit none

  integer fixed_num
  integer test

  if ( test == 1 ) then
    call p01_fixed_num ( fixed_num )
  else if ( test == 2 ) then
    call p02_fixed_num ( fixed_num )
  else if ( test == 3 ) then
    call p03_fixed_num ( fixed_num )
  else if ( test == 4 ) then
    call p04_fixed_num ( fixed_num )
  else if ( test == 5 ) then
    call p05_fixed_num ( fixed_num )
  else if ( test == 6 ) then
    call p06_fixed_num ( fixed_num )
  else if ( test == 7 ) then
    call p07_fixed_num ( fixed_num )
  else if ( test == 8 ) then
    call p08_fixed_num ( fixed_num )
  else if ( test == 9 ) then
    call p09_fixed_num ( fixed_num )
  else if ( test == 10 ) then
    call p10_fixed_num ( fixed_num )
  else if ( test == 11 ) then
    call p11_fixed_num ( fixed_num )
  else if ( test == 12 ) then
    call p12_fixed_num ( fixed_num )
  else if ( test == 13 ) then
    call p13_fixed_num ( fixed_num )
  else if ( test == 14 ) then
    call p14_fixed_num ( fixed_num )
  else if ( test == 15 ) then
    call p15_fixed_num ( fixed_num )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_FIXED_NUM - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_fixed_points ( test, m, fixed_num, fixed )

!*****************************************************************************80
!
!! P00_FIXED_POINTS returns the fixed points in a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer FIXED_NUM, the number of fixed points.
!
!    Output, real ( kind = 8 ) FIXED(M,FIXED_NUM), the 
!    coordinates of the fixed points.
!
  implicit none

  integer m
  integer fixed_num

  real ( kind = 8 ) fixed(m,fixed_num)
  integer test

  if ( test == 1 ) then
    call p01_fixed_points ( m, fixed_num, fixed )
  else if ( test == 2 ) then
    call p02_fixed_points ( m, fixed_num, fixed )
  else if ( test == 3 ) then
    call p03_fixed_points ( m, fixed_num, fixed )
  else if ( test == 4 ) then
    call p04_fixed_points ( m, fixed_num, fixed )
  else if ( test == 5 ) then
    call p05_fixed_points ( m, fixed_num, fixed )
  else if ( test == 6 ) then
    call p06_fixed_points ( m, fixed_num, fixed )
  else if ( test == 7 ) then
    call p07_fixed_points ( m, fixed_num, fixed )
  else if ( test == 8 ) then
    call p08_fixed_points ( m, fixed_num, fixed )
  else if ( test == 9 ) then
    call p09_fixed_points ( m, fixed_num, fixed )
  else if ( test == 10 ) then
    call p10_fixed_points ( m, fixed_num, fixed )
  else if ( test == 11 ) then
    call p11_fixed_points ( m, fixed_num, fixed )
  else if ( test == 12 ) then
    call p12_fixed_points ( m, fixed_num, fixed )
  else if ( test == 13 ) then
    call p13_fixed_points ( m, fixed_num, fixed )
  else if ( test == 14 ) then
    call p14_fixed_points ( m, fixed_num, fixed )
  else if ( test == 15 ) then
    call p15_fixed_points ( m, fixed_num, fixed )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_FIXED_POINTS - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_header ( test )

!*****************************************************************************80
!
!! P00_HEADER prints some information about a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
  implicit none

  integer test

  if ( test == 1 ) then
    call p01_header ( )
  else if ( test == 2 ) then
    call p02_header ( )
  else if ( test == 3 ) then
    call p03_header ( )
  else if ( test == 4 ) then
    call p04_header ( )
  else if ( test == 5 ) then
    call p05_header ( )
  else if ( test == 6 ) then
    call p06_header ( )
  else if ( test == 7 ) then
    call p07_header ( )
  else if ( test == 8 ) then
    call p08_header ( )
  else if ( test == 9 ) then
    call p09_header ( )
  else if ( test == 10 ) then
    call p10_header ( )
  else if ( test == 11 ) then
    call p11_header ( )
  else if ( test == 12 ) then
    call p12_header ( )
  else if ( test == 13 ) then
    call p13_header ( )
  else if ( test == 14 ) then
    call p14_header ( )
  else if ( test == 15 ) then
    call p15_header ( )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_HEADER - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_hex_grid ( test, m, h, n, point )

!*****************************************************************************80
!
!! P00_HEX_GRID returns hex grid points in a region.
!
!  Modified:
!
!    15 May 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer M, the spatial dimension.
!
!    Input, real ( kind = 8 ) H, the hexagonal spacing.
!
!    Input, integer N, the number of hexagonal grid points
!    that lie within the region, as computed by P00_HEX_GRID_COUNT.
!
!    Output, real POINT(M,N), the hex grid points.
!
  implicit none

  integer m

  real ( kind = 8 ) angle
  real ( kind = 8 ) box(m,2)
  real ( kind = 8 ) center(m)
  real ( kind = 8 ) h
  real ( kind = 8 ) hi(m)
  logical, allocatable, dimension ( : ) :: inside
  integer j
  integer j2
  real ( kind = 8 ) lo(m)
  integer n
  integer n2
  real ( kind = 8 ) point(m,n)
  real ( kind = 8 ), allocatable, dimension ( :, : ) :: point2
  integer test
!
!  Get the box limits.
!
  call p00_box ( test, m, lo, hi )
!
!  How many hex points will fit in the box?
!
  box(1:2,1) = lo(1:2)
  box(1:2,2) = hi(1:2)
  center(1:2) = 0.5D+00 * ( lo(1:2) + hi(1:2) )
  angle = 0.0D+00

  call hex_grid_angle_size ( box, center, angle, h, n2 )
!
!  Generate the hex points.
!
  allocate ( inside(1:n2) )
  allocate ( point2(1:m,1:n2) )

  call hex_grid_angle ( box, center, angle, h, n2, point2 )
!
!  How many of these points are in the region?
!
  call p00_inside ( test, m, n2, point2, inside )
!
!  Copy the good hex grid points.
!
  j = 0
  do j2 = 1, n2
    if ( inside(j2) ) then
      j = j + 1
      point(1:m,j) = point2(1:m,j2)
    end if
  end do

  deallocate ( inside )
  deallocate ( point2 )

  return
end
subroutine p00_hex_grid_count ( test, m, h, n )

!*****************************************************************************80
!
!! P00_HEX_GRID_COUNT counts the number of hex grid points in a region.
!
!  Modified:
!
!    15 May 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer M, the spatial dimension.
!
!    Input, real ( kind = 8 ) H, the hexagonal spacing.
!
!    Output, integer N, the number of hexagonal grid points
!    that lie within the region.
!
  implicit none

  integer m

  real ( kind = 8 ) angle
  real ( kind = 8 ) box(m,2)
  real ( kind = 8 ) center(m)
  real ( kind = 8 ) h
  real ( kind = 8 ) hi(m)
  logical, allocatable, dimension ( : ) :: inside
  integer j
  real ( kind = 8 ) lo(m)
  integer n
  integer n2
  real ( kind = 8 ), allocatable, dimension ( :, : ) :: point2
  integer test
!
!  Get the box limits.
!
  call p00_box ( test, m, lo, hi )
!
!  How many hex points will fit in the box?
!
  box(1:2,1) = lo(1:2)
  box(1:2,2) = hi(1:2)
  center(1:2) = 0.5D+00 * ( lo(1:2) + hi(1:2) )
  angle = 0.0D+00

  call hex_grid_angle_size ( box, center, angle, h, n2 )
!
!  Generate the hex points.
!
  allocate ( inside(1:n2) )
  allocate ( point2(1:m,1:n2) )

  call hex_grid_angle ( box, center, angle, h, n2, point2 )
!
!  How many of these points are in the region?
!
  call p00_inside ( test, m, n2, point2, inside )
!
!  Add them up.
!
  n = 0
  do j = 1, n2
    if ( inside(j) ) then
      n = n + 1
    end if
  end do

  deallocate ( inside )
  deallocate ( point2 )

  return
end
subroutine p00_hole_num ( test, hole_num )

!*****************************************************************************80
!
!! P00_HOLE_NUM counts the holes in a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Output, integer HOLE_NUM, the number of boundary segments.
!
  implicit none

  integer hole_num
  integer test

  if ( test == 1 ) then
    call p01_hole_num ( hole_num )
  else if ( test == 2 ) then
    call p02_hole_num ( hole_num )
  else if ( test == 3 ) then
    call p03_hole_num ( hole_num )
  else if ( test == 4 ) then
    call p04_hole_num ( hole_num )
  else if ( test == 5 ) then
    call p05_hole_num ( hole_num )
  else if ( test == 6 ) then
    call p06_hole_num ( hole_num )
  else if ( test == 7 ) then
    call p07_hole_num ( hole_num )
  else if ( test == 8 ) then
    call p08_hole_num ( hole_num )
  else if ( test == 9 ) then
    call p09_hole_num ( hole_num )
  else if ( test == 10 ) then
    call p10_hole_num ( hole_num )
  else if ( test == 11 ) then
    call p11_hole_num ( hole_num )
  else if ( test == 12 ) then
    call p12_hole_num ( hole_num )
  else if ( test == 13 ) then
    call p13_hole_num ( hole_num )
  else if ( test == 14 ) then
    call p14_hole_num ( hole_num )
  else if ( test == 15 ) then
    call p15_hole_num ( hole_num )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_HOLE_NUM - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    write ( *, '(a,i8)' ) '  TEST = ', test
    stop
  end if

  return
end
subroutine p00_hole_point ( test, hole_index, m, hole_point )

!*****************************************************************************80
!
!! P00_HOLE_POINT returns a point inside a given hole.
!
!  Modified:
!
!    01 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer HOLE_INDEX, the index of the hole.
!
!    Input, integer M, the spatial dimension.
!
!    Output, real ( kind = 8 ) HOLE_POINT(M), a point in the hole
!
  implicit none

  integer m

  integer hole_index
  real ( kind = 8 ) hole_point(m)
  integer test

  if ( test == 1 ) then
    call p01_hole_point ( hole_index, m, hole_point )
  else if ( test == 2 ) then
    call p02_hole_point ( hole_index, m, hole_point )
  else if ( test == 3 ) then
    call p03_hole_point ( hole_index, m, hole_point )
  else if ( test == 4 ) then
    call p04_hole_point ( hole_index, m, hole_point )
  else if ( test == 5 ) then
    call p05_hole_point ( hole_index, m, hole_point )
  else if ( test == 6 ) then
    call p06_hole_point ( hole_index, m, hole_point )
  else if ( test == 7 ) then
    call p07_hole_point ( hole_index, m, hole_point )
  else if ( test == 8 ) then
    call p08_hole_point ( hole_index, m, hole_point )
  else if ( test == 9 ) then
    call p09_hole_point ( hole_index, m, hole_point )
  else if ( test == 10 ) then
    call p10_hole_point ( hole_index, m, hole_point )
  else if ( test == 11 ) then
    call p11_hole_point ( hole_index, m, hole_point )
  else if ( test == 12 ) then
    call p12_hole_point ( hole_index, m, hole_point )
  else if ( test == 13 ) then
    call p13_hole_point ( hole_index, m, hole_point )
  else if ( test == 14 ) then
    call p14_hole_point ( hole_index, m, hole_point )
  else if ( test == 15 ) then
    call p15_hole_point ( hole_index, m, hole_point )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_HOLE_POINT - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_inside ( test, m, n, point, inside )

!*****************************************************************************80
!
!! P00_INSIDE reports if a point is inside the region in a problem.
!
!  Discussion:
!
!    For argument's sake, a point on the boundary can be considered
!    inside the region, but it is probably futile to attempt to distinguish
!    this case in general.  For more information about a point's location,
!    use P00_SDIST.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
!    Output, logical INSIDE(N), is TRUE if the point is in the region.
!
  implicit none

  integer m
  integer n

  logical inside(n)
  real ( kind = 8 ) point(m,n)
  integer test

  if ( test == 1 ) then
    call p01_inside ( m, n, point, inside )
  else if ( test == 2 ) then
    call p02_inside ( m, n, point, inside )
  else if ( test == 3 ) then
    call p03_inside ( m, n, point, inside )
  else if ( test == 4 ) then
    call p04_inside ( m, n, point, inside )
  else if ( test == 5 ) then
    call p05_inside ( m, n, point, inside )
  else if ( test == 6 ) then
    call p06_inside ( m, n, point, inside )
  else if ( test == 7 ) then
    call p07_inside ( m, n, point, inside )
  else if ( test == 8 ) then
    call p08_inside ( m, n, point, inside )
  else if ( test == 9 ) then
    call p09_inside ( m, n, point, inside )
  else if ( test == 10 ) then
    call p10_inside ( m, n, point, inside )
  else if ( test == 11 ) then
    call p11_inside ( m, n, point, inside )
  else if ( test == 12 ) then
    call p12_inside ( m, n, point, inside )
  else if ( test == 13 ) then
    call p13_inside ( m, n, point, inside )
  else if ( test == 14 ) then
    call p14_inside ( m, n, point, inside )
  else if ( test == 15 ) then
    call p15_inside ( m, n, point, inside )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_INSIDE - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_points_eps ( test, h, m, n, points, eps_file_name )

!*****************************************************************************80
!
!! P00_POINTS_EPS draws points in a region as an EPS file.
!
!  Discussion:
!
!    The boundary of the region is also drawn.
!
!  Modified:
!
!    09 June 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the number of the problem.
!
!    Input, real ( kind = 8 ) H, the maximum size of a segment of the boundary.
!    This controls how smoothly curved sections of the boundary will be drawn.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, integer POINTS(M,N), the points to draw.
!
!    Input, character ( len = * ) EPS_FILE_NAME, the name of the EPS file
!    to create.
!
!  Local Parameters:
!
!    Local, integer CIRCLE_SIZE, controls the size of the circles
!    used to indicate points.  These are measured in PostScript points,
!    that is, 1/72 of an inch.  A value of 3, 4 or 5 is usually reasonable.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ), allocatable, dimension ( :, : ) :: boundary
  integer, parameter :: circle_size = 4
  character ( len = * ) eps_file_name
  integer eps_unit
  real ( kind = 8 ) h
  real ( kind = 8 ), allocatable, dimension ( : ) :: hi
  integer ios
  integer j
  real ( kind = 8 ), allocatable, dimension ( : ) :: lo
  real ( kind = 8 ) points(m,n)
  real ( kind = 8 ) scale
  integer segment
  integer segment_length
  integer segment_num
  logical, parameter :: show_points = .false.
  character ( len = 40 ) string
  integer test
  real ( kind = 8 ) x_max
  real ( kind = 8 ) x_min
  integer x_ps
  integer :: x_ps_max = 576
  integer :: x_ps_max_clip = 594
  integer :: x_ps_max_user
  integer :: x_ps_min = 36
  integer :: x_ps_min_clip = 18
  integer :: x_ps_min_user
  real ( kind = 8 ) x_scale
  real ( kind = 8 ) y_max
  real ( kind = 8 ) y_min
  integer y_ps
  integer :: y_ps_max = 666
  integer :: y_ps_max_clip = 684
  integer :: y_ps_max_user
  integer :: y_ps_min = 126
  integer :: y_ps_min_clip = 108
  integer :: y_ps_min_user
  real ( kind = 8 ) y_scale

  call p00_boundary_segment_num ( test, segment_num )

  allocate ( lo(1:m) )
  allocate ( hi(1:m) )

  call p00_box ( test, m, lo, hi )

  x_min = lo(1) - 0.025D+00 * ( hi(1) - lo(1) )
  y_min = lo(2) - 0.025D+00 * ( hi(2) - lo(2) )
  x_max = hi(1) + 0.025D+00 * ( hi(1) - lo(1) )
  y_max = hi(2) + 0.025D+00 * ( hi(2) - lo(2) )

  x_min = min ( x_min, minval ( points(1,1:n) ) )
  x_max = max ( x_max, maxval ( points(1,1:n) ) )
  y_min = min ( y_min, minval ( points(2,1:n) ) )
  y_max = max ( y_max, maxval ( points(2,1:n) ) )

  x_scale = x_max - x_min
  y_scale = y_max - y_min
  scale = max ( x_scale, y_scale )
!
!  Determine the PostScript coordinates of the used box.
!
  x_ps_min_user = ( ( x_ps_max + x_ps_min ) &
    - int ( real ( x_ps_max - x_ps_min, kind = 8 ) * x_scale / scale ) ) / 2
  x_ps_max_user = ( ( x_ps_max + x_ps_min ) &
    + int ( real ( x_ps_max - x_ps_min, kind = 8 ) * x_scale / scale ) ) / 2
  y_ps_min_user = ( ( y_ps_max + y_ps_min ) &
    - int ( real ( y_ps_max - y_ps_min, kind = 8 ) * y_scale / scale ) ) / 2
  y_ps_max_user = ( ( y_ps_max + y_ps_min ) &
    + int ( real ( y_ps_max - y_ps_min, kind = 8 ) * y_scale / scale ) ) / 2

  if ( x_scale < y_scale ) then
    x_max = x_max + 0.5D+00 * ( y_scale - x_scale )
    x_min = x_min - 0.5D+00 * ( y_scale - x_scale )
  else if ( y_scale < x_scale ) then
    y_max = y_max + 0.5D+00 * ( x_scale - y_scale )
    y_min = y_min - 0.5D+00 * ( x_scale - y_scale )
  end if

  call get_unit ( eps_unit )

  open ( unit = eps_unit, file = eps_file_name, status = 'replace', &
    iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_POINTS_EPS - Fatal error!'
    write ( *, '(a)' ) '  Could not open the output EPS file.'
    stop
  end if

  call timestring ( string )

  write ( eps_unit, '(a)' ) '%!PS-Adobe-3.0 EPSF-3.0'
  write ( eps_unit, '(a)' ) '%%Creator: triangulation_plot_eps.f90'
  write ( eps_unit, '(a)' ) '%%Title: ' // trim ( eps_file_name )
  write ( eps_unit, '(a)' ) '%%CreationDate: ' // trim ( string )
  write ( eps_unit, '(a)' ) '%%Pages: 1'
  write ( eps_unit, '(a,i4,2x,i4,2x,i4,2x,i4)' ) '%%BoundingBox: ', &
    x_ps_min_user, y_ps_min_user, x_ps_max_user, y_ps_max_user
  write ( eps_unit, '(a)' ) '%%Document-Fonts: Times-Roman'
  write ( eps_unit, '(a)' ) '%%LanguageLevel: 1'
  write ( eps_unit, '(a)' ) '%%EndComments'
  write ( eps_unit, '(a)' ) '%%BeginProlog'
  write ( eps_unit, '(a)' ) '/inch {72 mul} def'
  write ( eps_unit, '(a)' ) '%%EndProlog'
  write ( eps_unit, '(a)' ) '%%Page: 1 1'
  write ( eps_unit, '(a)' ) 'save'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Set the RGB line color to very light gray.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '0.900  0.900  0.900 setrgbcolor'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Draw a gray border around the user box.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) 'newpath'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_min_user, y_ps_min_user, ' moveto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_max_user, y_ps_min_user, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_max_user, y_ps_max_user, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_min_user, y_ps_max_user, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) &
    '  ', x_ps_min_user, y_ps_min_user, ' lineto'
  write ( eps_unit, '(a)' ) 'stroke'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Set the RGB line color to black.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '0.000  0.000  0.000 setrgbcolor'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Set the font and its size.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '/Times-Roman findfont'
  write ( eps_unit, '(a)' ) '0.50 inch scalefont'
  write ( eps_unit, '(a)' ) 'setfont'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Define a clipping polygon.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) 'newpath'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_min_clip, y_ps_min_clip, ' moveto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_max_clip, y_ps_min_clip, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_max_clip, y_ps_max_clip, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_min_clip, y_ps_max_clip, ' lineto'
  write ( eps_unit, '(a,i4,2x,i4,2x,a)' ) '  ', &
    x_ps_min_clip, y_ps_min_clip, ' lineto'
  write ( eps_unit, '(a)' ) 'clip newpath'

  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Set the RGB line color to red.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '0.900  0.200  0.100 setrgbcolor'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Increase the linewidth to 3.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '3  setlinewidth'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Draw the boundary lines.'
  write ( eps_unit, '(a)' ) '%'

  do segment = 1, segment_num

    call p00_boundary_segment_length ( test, segment, h, segment_length )

    if ( segment_length <= 0 ) then
      cycle
    end if

    allocate ( boundary(1:m,1:segment_length) )

    call p00_boundary_segment ( test, segment, m, segment_length, boundary )

    do j = 1, segment_length

      x_ps = int ( &
        ( ( x_max - boundary(1,j)         ) * real ( x_ps_min, kind = 8 )   &
        + (         boundary(1,j) - x_min ) * real ( x_ps_max, kind = 8 ) ) &
        / ( x_max                 - x_min ) )

      y_ps = int ( &
        ( ( y_max - boundary(2,j)         ) * real ( y_ps_min, kind = 8 )   &
        + (         boundary(2,j) - y_min ) * real ( y_ps_max, kind = 8 ) ) &
        / ( y_max                 - y_min ) )

      if ( j == 1 ) then
        write ( eps_unit, '(i4,2x,i4,2x,a)' ) x_ps, y_ps, ' moveto'
      else
        write ( eps_unit, '(i4,2x,i4,2x,a)' ) x_ps, y_ps, ' lineto'
      end if

    end do

    write ( eps_unit, '(a)' ) 'stroke'

    deallocate ( boundary )

  end do

  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Set the RGB line color to blue.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '0.000  0.150  0.750 setrgbcolor'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Restore the linewidth to 1.'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '1  setlinewidth'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%  Draw the nodes.'
  write ( eps_unit, '(a)' ) '%'

  do j = 1, n

    x_ps = int ( &
      ( ( x_max - points(1,j)         ) * real ( x_ps_min, kind = 8 )   &
      + (         points(1,j) - x_min ) * real ( x_ps_max, kind = 8 ) ) &
      / ( x_max               - x_min ) )

    y_ps = int ( &
      ( ( y_max - points(2,j)         ) * real ( y_ps_min, kind = 8 )   &
      + (         points(2,j) - y_min ) * real ( y_ps_max, kind = 8 ) ) &
      / ( y_max               - y_min ) )

    write ( eps_unit, '(a,i4,2x,i4,2x,i4,2x,a)' ) 'newpath ', x_ps, y_ps, &
      circle_size, ' 0 360 arc closepath fill'

  end do

  deallocate ( hi )
  deallocate ( lo )

  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) 'restore showpage'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '% End of page'
  write ( eps_unit, '(a)' ) '%'
  write ( eps_unit, '(a)' ) '%%Trailer'
  write ( eps_unit, '(a)' ) '%%EOF'

  close ( unit = eps_unit )

  if ( .false. ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_POINTS_EPS:'
    write ( *, '(a)' ) '  An encapsulated PostScript file was created'
    write ( *, '(a)' ) '  containing an image of the points.'
    write ( *, '(a)' ) '  The file is named "' // trim ( eps_file_name ) // '".'
  end if

  return
end
subroutine p00_poly_write ( test, file_name )

!*****************************************************************************80
!
!! P00_POLY_WRITE collects data and writes it to a POLY file.
!
!  Discussion:
!
!    This routine collects information about the boundary for a given
!    problem, and writes that data to a POLY file that can be read by
!    TRIANGLE.
!
!  Modified:
!
!    03 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, character ( len = * ) FILE_NAME, the name of the file to create.
!
  implicit none

  integer, parameter :: dim_num = 2

  integer boundary_segment_num
  integer e
  integer, allocatable, dimension ( :, : ) :: edge_nodes
  integer edge_num
  character ( len = * ) file_name
  real ( kind = 8 ) h
  real ( kind = 8 ) hi(dim_num)
  integer hole_index
  integer hole_num
  real ( kind = 8 ), allocatable, dimension ( :, : ) :: hole
  integer j
  real ( kind = 8 ) lo(dim_num)
  integer next
  integer node_num
  real ( kind = 8 ) scale
  real ( kind = 8 ), allocatable, dimension ( :, : ) :: segment
  integer segment_index
  integer segment_length
  integer test
!
!  Get a length scale.
!
  call p00_box ( test, dim_num, lo, hi )

  scale = maxval ( abs ( hi(1:dim_num) - lo(1:dim_num) ) )
!
!  How many boundary segments are there?
!
  call p00_boundary_segment_num ( test, boundary_segment_num )
!
!  Choose H so that we would expect to get about 100 boundary points if our
!  region were a square of any size.
!
  h = 0.04D+00 * scale

  node_num = 0
  do segment_index = 1, boundary_segment_num
    call p00_boundary_segment_length ( test, segment_index, h, segment_length )
    node_num = node_num + segment_length
  end do

  allocate ( segment(dim_num,node_num) )
!
!  Now collect all the boundary nodes into one array.
!
  next = 1

  do segment_index = 1, boundary_segment_num

    call p00_boundary_segment_length ( test, segment_index, h, segment_length )

    call p00_boundary_segment ( test, segment_index, dim_num, &
      segment_length, segment(1:dim_num,next:next+segment_length-1) )

    next = next + segment_length

  end do
!
!  How many edges are there?
!
  edge_num = 0

  do segment_index = 1, boundary_segment_num

    call p00_boundary_segment_length ( test, segment_index, h, segment_length )

    edge_num = edge_num + segment_length - 1

  end do

  allocate ( edge_nodes(2,edge_num) )
!
!  Now collect the edges.
!
  e = 0
  next = 1

  do segment_index = 1, boundary_segment_num

    call p00_boundary_segment_length ( test, segment_index, h, segment_length )

    do j = 1, segment_length - 1
      edge_nodes(1,e+j) = next + j - 1
      edge_nodes(2,e+j) = next + j
    end do

    next = next + segment_length
    e = e + segment_length - 1

  end do
!
!  Handle the holes.
!
  call p00_hole_num ( test, hole_num )

  allocate ( hole(dim_num,hole_num) )

  do hole_index = 1, hole_num
    call p00_hole_point ( test, hole_index, dim_num, &
      hole(1:dim_num,hole_index) )
  end do
!
!  Write the POLY file.
!
  call poly_write ( file_name, node_num, segment, edge_num, &
    edge_nodes, hole_num, hole )

  deallocate ( edge_nodes )
  deallocate ( hole )
  deallocate ( segment )

  return
end
subroutine p00_sample ( test, m, n, seed, point )

!*****************************************************************************80
!
!! P00_SAMPLE samples points from the region in a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input/output, integer SEED, a seed for the random number generator.
!
!    Output, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ) point(m,n)
  integer seed
  integer test

  if ( test == 1 ) then
    call p01_sample ( m, n, seed, point )
  else if ( test == 2 ) then
    call p02_sample ( m, n, seed, point )
  else if ( test == 3 ) then
    call p03_sample ( m, n, seed, point )
  else if ( test == 4 ) then
    call p04_sample ( m, n, seed, point )
  else if ( test == 5 ) then
    call p05_sample ( m, n, seed, point )
  else if ( test == 6 ) then
    call p06_sample ( m, n, seed, point )
  else if ( test == 7 ) then
    call p07_sample ( m, n, seed, point )
  else if ( test == 8 ) then
    call p08_sample ( m, n, seed, point )
  else if ( test == 9 ) then
    call p09_sample ( m, n, seed, point )
  else if ( test == 10 ) then
    call p10_sample ( m, n, seed, point )
  else if ( test == 11 ) then
    call p11_sample ( m, n, seed, point )
  else if ( test == 12 ) then
    call p12_sample ( m, n, seed, point )
  else if ( test == 13 ) then
    call p13_sample ( m, n, seed, point )
  else if ( test == 14 ) then
    call p14_sample ( m, n, seed, point )
  else if ( test == 15 ) then
    call p15_sample ( m, n, seed, point )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_SAMPLE - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_sample_h1 ( test, m, n, h, seed, point )

!*****************************************************************************80
!
!! P00_SAMPLE_H1 samples points from the enlarged region in a problem.
!
!  Discussion:
!
!    The region is enlarged by an amount H.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) H, the enlargement amount.
!
!    Input/output, integer SEED, a seed for the random number generator.
!
!    Output, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ) h
  real ( kind = 8 ) point(m,n)
  integer seed
  integer test

  if ( test == 1 ) then
    call p01_sample_h1 ( m, n, h, seed, point )
  else if ( test == 2 ) then
    call p02_sample_h1 ( m, n, h, seed, point )
  else if ( test == 3 ) then
    call p03_sample_h1 ( m, n, h, seed, point )
  else if ( test == 4 ) then
    call p04_sample_h1 ( m, n, h, seed, point )
  else if ( test == 5 ) then
   call p05_sample_h1 ( m, n, h, seed, point )
  else if ( test == 6 ) then
    call p06_sample_h1 ( m, n, h, seed, point )
  else if ( test == 7 ) then
    call p07_sample_h1 ( m, n, h, seed, point )
  else if ( test == 8 ) then
    call p08_sample_h1 ( m, n, h, seed, point )
  else if ( test == 9 ) then
    call p09_sample_h1 ( m, n, h, seed, point )
  else if ( test == 10 ) then
    call p10_sample_h1 ( m, n, h, seed, point )
  else if ( test == 11 ) then
    call p11_sample_h1 ( m, n, h, seed, point )
  else if ( test == 12 ) then
    call p12_sample_h1 ( m, n, h, seed, point )
  else if ( test == 13 ) then
    call p13_sample_h1 ( m, n, h, seed, point )
  else if ( test == 14 ) then
    call p14_sample_h1 ( m, n, h, seed, point )
  else if ( test == 15 ) then
    call p15_sample_h1 ( m, n, h, seed, point )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_SAMPLE_H1 - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_sdist ( test, m, n, point, sdist )

!*****************************************************************************80
!
!! P00_SDIST returns the signed distance to the region in a problem.
!
!  Discussion:
!
!    A positive distance indicates the point is outside the region.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
!    Output, real ( kind = 8 ) SDIST(N), the signed distance 
!    of each point to the region.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ) point(m,n)
  real ( kind = 8 ) sdist(n)
  integer test

  if ( test == 1 ) then
    call p01_sdist ( m, n, point, sdist )
  else if ( test == 2 ) then
    call p02_sdist ( m, n, point, sdist )
  else if ( test == 3 ) then
    call p03_sdist ( m, n, point, sdist )
  else if ( test == 4 ) then
    call p04_sdist ( m, n, point, sdist )
  else if ( test == 5 ) then
    call p05_sdist ( m, n, point, sdist )
  else if ( test == 6 ) then
    call p06_sdist ( m, n, point, sdist )
  else if ( test == 7 ) then
    call p07_sdist ( m, n, point, sdist )
  else if ( test == 8 ) then
    call p08_sdist ( m, n, point, sdist )
  else if ( test == 9 ) then
    call p09_sdist ( m, n, point, sdist )
  else if ( test == 10 ) then
    call p10_sdist ( m, n, point, sdist )
  else if ( test == 11 ) then
    call p11_sdist ( m, n, point, sdist )
  else if ( test == 12 ) then
    call p12_sdist ( m, n, point, sdist )
  else if ( test == 13 ) then
    call p13_sdist ( m, n, point, sdist )
  else if ( test == 14 ) then
    call p14_sdist ( m, n, point, sdist )
  else if ( test == 15 ) then
    call p15_sdist ( m, n, point, sdist )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_SDIST - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p00_test_num ( test_num )

!*****************************************************************************80
!
!! P00_TEST_NUM returns the number of available tests.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer TEST_NUM, the number of tests.
!
  implicit none

  integer test_num

  test_num = 15

  return
end
subroutine p00_title ( test, title )

!*****************************************************************************80
!
!! P00_TITLE returns a title for a problem.
!
!  Modified:
!
!    06 February 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer TEST, the index of the test problem
!
!    Output, character ( len = * ) TITLE, a title for the problem.
!
  implicit none

  integer test
  character ( len = * ) title

  if ( test == 1 ) then
    call p01_title ( title )
  else if ( test == 2 ) then
    call p02_title ( title )
  else if ( test == 3 ) then
    call p03_title ( title )
  else if ( test == 4 ) then
    call p04_title ( title )
  else if ( test == 5 ) then
    call p05_title ( title )
  else if ( test == 6 ) then
    call p06_title ( title )
  else if ( test == 7 ) then
    call p07_title ( title )
  else if ( test == 8 ) then
    call p08_title ( title )
  else if ( test == 9 ) then
    call p09_title ( title )
  else if ( test == 10 ) then
    call p10_title ( title )
  else if ( test == 11 ) then
    call p11_title ( title )
  else if ( test == 12 ) then
    call p12_title ( title )
  else if ( test == 13 ) then
    call p13_title ( title )
  else if ( test == 14 ) then
    call p14_title ( title )
  else if ( test == 15 ) then
    call p15_title ( title )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_TITLE - Fatal error!'
    write ( *, '(a)' ) '  Input value of TEST is out of range.'
    stop
  end if

  return
end
subroutine p01_boundary_nearest ( m, n, point, boundary )

!*****************************************************************************80
!
!! P01_BOUNDARY_NEAREST returns a nearest boundary point in problem 01.
!
!  Discussion:
!
!    The given input point need not be inside the region.
!
!    In some cases, more than one boundary point may be "nearest",
!    but only one will be returned.
!
!  Modified:
!
!    01 October 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
!    Output, real ( kind = 8 ) BOUNDARY(M,N), points on the boundary
!    that are nearest to each point.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ), dimension ( m, n ) :: boundary
  real ( kind = 8 ), dimension ( 2 ) :: center =  (/ &
    0.0D+00, 0.0D+00 /)
  integer j
  real ( kind = 8 ), dimension ( m, n ) :: point
  real ( kind = 8 ) :: r
  real ( kind = 8 ), parameter :: r1 = 1.0D+00

  do j = 1, n

    r = sqrt ( sum ( ( point(1:m,j) - center(1:m) )**2 ) )

    if ( r == 0.0D+00 ) then
      boundary(1,j) = center(1) + r1
      boundary(2,j) = center(2)
    else
      boundary(1:m,j) = center(1:m) &
        + ( r1 / r ) * ( point(1:m,j) - center(1:m) )
    end if

  end do

  return
end
subroutine p01_boundary_project ( m, n, point )

!*****************************************************************************80
!
!! P01_BOUNDARY_PROJECT projects exterior points to the boundary in problem 01.
!
!  Discussion:
!
!    Points in the interior are not changed.
!
!  Modified:
!
!    11 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input/output, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.  On output, all exterior points have been
!    replaced by the nearest point on the boundary.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ), dimension ( 2 ), parameter :: center =  &
    (/ 0.0D+00, 0.0D+00 /)
  logical inside(n)
  integer j
  real ( kind = 8 ), dimension ( m, n ) :: point
  real ( kind = 8 ) :: r
  real ( kind = 8 ), parameter :: r1 = 1.0D+00

  call p01_inside ( m, n, point, inside )

  do j = 1, n

    if ( inside(j) ) then
      cycle
    end if

    r = sqrt ( sum ( ( point(1:m,j) - center(1:m) )**2 ) )

    if ( r == 0.0D+00 ) then
      point(1,j) = center(1) + r1
      point(2,j) = center(2)
    else
      point(1:m,j) = center(1:m) &
        + ( r1 / r ) * ( point(1:m,j) - center(1:m) )
    end if

  end do

  return
end
subroutine p01_boundary_segment ( segment_index, m, segment_length, &
  segment )

!*****************************************************************************80
!
!! P01_BOUNDARY_SEGMENT returns a boundary segment in problem 01.
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer SEGMENT_INDEX, the index of the boundary segment.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer SEGMENT_LENGTH, the number of points in the segment.
!
!    Output, real ( kind = 8 ) SEGMENT(M,SEGMENT_LENGTH), the 
!    points that make up the boundary segment.
!
  implicit none

  integer m
  integer segment_length

  real ( kind = 8 ) angle
  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  integer i
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  integer segment_index
  real ( kind = 8 ) segment(m,segment_length)

  if ( segment_index == 1 ) then

    do i = 1, segment_length

      angle = 2.0D+00 * pi &
        * real (              i - 1, kind = 8 ) &
        / real ( segment_length - 1, kind = 8 )

      segment(1,i) = center(1) + r1 * cos ( angle )
      segment(2,i) = center(2) + r1 * sin ( angle )

    end do

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P01_BOUNDARY_SEGMENT - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal SEGMENT_INDEX = ', segment_index
    stop

  end if

  return
end
subroutine p01_boundary_segment_length ( segment_index, h, segment_length )

!*****************************************************************************80
!
!! P01_BOUNDARY_SEGMENT_LENGTH returns boundary segment lengths in problem 01.
!
!  Modified:
!
!    13 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer SEGMENT_INDEX, the index of one of the boundary segments.
!
!    Input, real ( kind = 8 ) H, the suggested spacing between points.
!
!    Output, integer SEGMENT_LENGTH, the number of points in the segment.
!
  implicit none

  real ( kind = 8 ) h
  integer n
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  integer segment_index
  integer segment_length

  if ( h <= 0.0D+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P01_BOUNDARY_SEGMENT_LENGTH - Fatal error!'
    write ( *, '(a,g14.6)' ) '  Nonpositive H = ', h
    stop
  end if

  if ( segment_index == 1 ) then

    n = nint ( 2.0D+00 * pi * r1 / h )
    n = max ( n, 5 )
    segment_length = n

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P01_BOUNDARY_SEGMENT_LENGTH - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal SEGMENT_INDEX = ', segment_index
    stop

  end if

  return
end
subroutine p01_boundary_segment_num ( boundary_segment_num )

!*****************************************************************************80
!
!! P01_BOUNDARY_SEGMENT_NUM counts the boundary segments in problem 01.
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Output, integer BOUNDARY_SEGMENT_NUM, the number of boundary segments.
!
  implicit none

  integer boundary_segment_num

  boundary_segment_num = 1

  return
end
subroutine p01_box ( m, lo, hi )

!*****************************************************************************80
!
!! P01_BOX returns a bounding box for problem 01.
!
!  Modified:
!
!    11 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Output, real ( kind = 8 ) LO(M), HI(M), coordinates of the
!    low and high corners of the box.
!
  implicit none

  integer m

  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  real ( kind = 8 ) hi(m)
  real ( kind = 8 ) lo(m)
  real ( kind = 8 ), parameter :: r1 = 1.0D+00

  lo(1:m) = (/ center(1) - r1, center(2) - r1 /)
  hi(1:m) = (/ center(1) + r1, center(2) + r1 /)

  return
end
subroutine p01_density ( m, n, point, density )

!*****************************************************************************80
!
!! P01_DENSITY returns the density for problem 01.
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
!    Output, real ( kind = 8 ) DENSITY(N), the mesh density 
!    at each point.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ) density(n)
  real ( kind = 8 ) point(m,n)

  density(1:n) = 1.0D+00

  return
end
subroutine p01_element_size ( element_size )

!*****************************************************************************80
!
!! P01_ELEMENT_SIZE returns a typical element size for problem 01.
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ELEMENT_SIZE, a typical element size.
!
  implicit none

  real ( kind = 8 ) element_size

  element_size = 0.2D+00

  return
end
subroutine p01_fixed_num ( fixed_num )

!*****************************************************************************80
!
!! P01_FIXED_NUM returns the number of fixed points in problem 01.
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Output, integer FIXED_NUM, the number of fixed points.
!
  implicit none

  integer fixed_num

  fixed_num = 0

  return
end
subroutine p01_fixed_points ( m, fixed_num, fixed )

!*****************************************************************************80
!
!! P01_FIXED_POINTS returns the fixed points in problem 01.
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer FIXED_NUM, the number of fixed points.
!
!    Output, real ( kind = 8 ) FIXED(M,FIXED_NUM), the coordinates 
!    of the fixed points.
!
  implicit none

  integer m
  integer fixed_num

  real ( kind = 8 ) fixed(m,fixed_num)

  return
end
subroutine p01_header ( )

!*****************************************************************************80
!
!! P01_HEADER prints some information about problem 01.
!
!  Modified:
!
!    11 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    None.
!
  implicit none

  integer, parameter :: m = 2

  integer boundary_segment_num
  real ( kind = 8 ), dimension ( m ) :: center = (/ 0.0D+00, 0.0D+00 /)
  integer fixed_num
  integer hole_num
  real ( kind = 8 ), parameter :: r1 = 1.0D+00

  call p01_boundary_segment_num ( boundary_segment_num ) 
  call p01_fixed_num ( fixed_num )  
  call p01_hole_num ( hole_num )  

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'P01:'
  write ( *, '(a)' ) '  Strang and Persson example #1'
  write ( *, '(a)' ) '  The unit circle.'
  write ( *, '(a,g14.6)' ) '  Radius = ', r1
  write ( *, '(a,2g14.6)' ) '  Center = ', center(1:2)
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  A uniform mesh density is requested.'
  write ( *, '(a)' ) '  Element sizes tried were 0.4, 0.2, 0.1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of boundary segments = ', boundary_segment_num
  write ( *, '(a,i8)' ) '  Number of fixed points =      ', fixed_num
  write ( *, '(a,i8)' ) '  Number of holes =             ', hole_num

  return
end
subroutine p01_hole_num ( hole_num )

!*****************************************************************************80
!
!! P01_HOLE_NUM counts the holes in problem 01.
!
!  Modified:
!
!    11 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer HOLE_NUM, the number of holes.
!
  implicit none

  integer hole_num

  hole_num = 0

  return
end
subroutine p01_hole_point ( hole_index, m, hole_point )

!*****************************************************************************80
!
!! P01_HOLE_POINT returns a point inside a given hole in problem 1.
!
!  Modified:
!
!    11 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer HOLE_INDEX, the index of the hole.
!
!    Input, integer M, the spatial dimension.
!
!    Output, real ( kind = 8 ) HOLE_POINT(M), a point in the hole
!
  implicit none

  integer m

  integer hole_index
  real ( kind = 8 ) hole_point(m)

  return
end
subroutine p01_inside ( m, n, point, inside )

!*****************************************************************************80
!
!! P01_INSIDE reports if a point is inside the region in problem 01.
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(M,N), the coordinates
!    of the points.
!
!    Output, logical INSIDE(N), is TRUE if the point is in the region.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  logical inside(n)
  real ( kind = 8 ) point(m,n)
  real ( kind = 8 ), parameter :: r1 = 1.0D+00

  inside(1:n) = ( ( point(1,1:n) - center(1) )**2 &
                + ( point(2,1:n) - center(2) )**2 ) <= r1**2

  return
end
subroutine p01_sample ( m, n, seed, point )

!*****************************************************************************80
!
!! P01_SAMPLE samples points from the region in problem 01.
!
!  Modified:
!
!    29 October 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input/output, integer SEED, a seed for the random number generator.
!
!    Output, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ) angle(n)
  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) point(m,n)
  real ( kind = 8 ) r(n)
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  integer seed
!
!  Choose uniform random angles between 0 and 2*Pi.
!
  call r8vec_uniform_01 ( n, seed, angle )

  angle(1:n) = 2.0D+00 * pi * angle(1:n)
!
!  Choose uniform random radii, then take square root.
!
  call r8vec_uniform_01 ( n, seed, r )

  r(1:n) = sqrt ( r(1:n) )
!
!  Construct the uniformly random points in the circle of radius R1 
!  centered at CENTER.
!
  point(1,1:n) = center(1) + r1 * r(1:n) * cos ( angle(1:n) )
  point(2,1:n) = center(2) + r1 * r(1:n) * sin ( angle(1:n) )

  return
end
subroutine p01_sample_h1 ( m, n, h, seed, point )

!*****************************************************************************80
!
!! P01_SAMPLE_H1 samples points from the enlarged region in problem 01.
!
!  Modified:
!
!    11 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) H, the enlargement amount.
!
!    Input/output, integer SEED, a seed for the random number generator.
!
!    Output, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ) angle(n)
  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  real ( kind = 8 ) h
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) point(m,n)
  real ( kind = 8 ) r(n)
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  integer seed
!
!  Choose uniform random angles between 0 and 2*Pi.
!
  call r8vec_uniform_01 ( n, seed, angle )

  angle(1:n) = 2.0D+00 * pi * angle(1:n)
!
!  Choose uniform random radii, then take square root.
!
  call r8vec_uniform_01 ( n, seed, r )

  r(1:n) = sqrt ( r(1:n) )
!
!  Construct the uniformly random points in the circle of radius R1 + H
!  centered at CENTER.
!
  point(1,1:n) = center(1) + ( r1 + h ) * r(1:n) * cos ( angle(1:n) )
  point(2,1:n) = center(2) + ( r1 + h ) * r(1:n) * sin ( angle(1:n) )

  return
end
subroutine p01_sdist ( m, n, point, sdist )

!*****************************************************************************80
!
!! P01_SDIST returns the signed distance to the region in problem 01.
!
!  Discussion:
!
!    A positive distance indicates the point is outside the region.
!
!  Modified:
!
!    11 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
!    Output, real ( kind = 8 ) SDIST(N), the signed distance of 
!    each point to the region.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  real ( kind = 8 ) point(m,n)
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  real ( kind = 8 ) sdist(n)

  sdist(1:n) = sqrt ( ( point(1,1:n) - center(1) )**2 &
                    + ( point(2,1:n) - center(2) )**2 ) - r1

  return
end
subroutine p01_title ( title )

!*****************************************************************************80
!
!! P01_TITLE returns a title for problem 01.
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Output, character ( len = * ) TITLE, a title for the problem.
!
  implicit none

  character ( len = * ) title

  title = '#1: The unit circle.'

  return
end
subroutine p02_boundary_nearest ( m, n, point, boundary )

!*****************************************************************************80
!
!! P02_BOUNDARY_NEAREST returns a nearest boundary point in problem 02.
!
!  Discussion:
!
!    The given input point need not be inside the region.
!
!    In some cases, more than one boundary point may be "nearest",
!    but only one will be returned.
!
!  Modified:
!
!    11 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
!    Output, real ( kind = 8 ) BOUNDARY(M,N), points on the boundary
!    that are nearest to each point.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ), dimension ( m, n ) :: boundary
  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  integer j
  real ( kind = 8 ), dimension ( m, n ) :: point
  real ( kind = 8 ) r
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  real ( kind = 8 ), parameter :: r2 = 0.4D+00

  do j = 1, n

    if ( point(1,j) == center(1) .and. point(2,j) == center(2) ) then

      boundary(1,j) = r2
      boundary(2,j) = center(2)

    else

      r = sqrt ( ( point(1,j) - center(1) )**2 &
               + ( point(2,j) - center(2) )**2 )

      if ( r1 - r < r - r2 ) then
        boundary(1,j) = center(1) + r1 * ( point(1,j) - center(1) ) / r
        boundary(2,j) = center(2) + r1 * ( point(2,j) - center(2) ) / r
      else
        boundary(1,j) = center(1) + r2 * ( point(1,j) - center(1) ) / r
        boundary(2,j) = center(2) + r2 * ( point(2,j) - center(2) ) / r
      end if

    end if

  end do

  return
end
subroutine p02_boundary_project ( m, n, point )

!*****************************************************************************80
!
!! P02_BOUNDARY_PROJECT projects exterior points to the boundary in problem 02.
!
!  Discussion:
!
!    Points in the interior are not changed.
!
!  Modified:
!
!    11 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input/output, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.  On output, all exterior points have been
!    replaced by the nearest point on the boundary.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  logical inside(n)
  integer j
  real ( kind = 8 ), dimension ( m, n ) :: point
  real ( kind = 8 ) r
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  real ( kind = 8 ), parameter :: r2 = 0.4D+00

  call p02_inside ( m, n, point, inside )

  do j = 1, n

    if ( inside(j) ) then
      cycle
    end if

    if ( point(1,j) == center(1) .and. point(2,j) == center(2) ) then

      point(1,j) = r2
      point(2,j) = center(2)

    else

      r = sqrt ( ( point(1,j) - center(1) )**2 &
               + ( point(2,j) - center(2) )**2 )

      if ( r1 - r < r - r2 ) then
        point(1,j) = center(1) + r1 * ( point(1,j) - center(1) ) / r
        point(2,j) = center(2) + r1 * ( point(2,j) - center(2) ) / r
      else
        point(1,j) = center(1) + r2 * ( point(1,j) - center(1) ) / r
        point(2,j) = center(2) + r2 * ( point(2,j) - center(2) ) / r
      end if

    end if

  end do

  return
end
subroutine p02_boundary_segment ( segment_index, m, segment_length, &
  segment )

!*****************************************************************************80
!
!! P02_BOUNDARY_SEGMENT returns a boundary segment in problem 02.
!
!  Modified:
!
!    15 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer SEGMENT_INDEX, the index of the boundary segment.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer SEGMENT_LENGTH, the number of points in the segment.
!
!    Output, real ( kind = 8 ) SEGMENT(M,SEGMENT_LENGTH), the 
!    points that make up the boundary segment.
!
  implicit none

  integer m
  integer segment_length

  real ( kind = 8 ) angle
  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  integer j
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  real ( kind = 8 ), parameter :: r2 = 0.4D+00
  integer segment_index
  real ( kind = 8 ) segment(m,segment_length)

  if ( segment_index == 1 ) then

    do j = 1, segment_length
      angle = 2.0D+00 * pi &
        * real (              j - 1, kind = 8 ) &
        / real ( segment_length - 1, kind = 8 )
      segment(1,j) = center(1) + r1 * cos ( angle )
      segment(2,j) = center(2) + r1 * sin ( angle )
    end do

  else if ( segment_index == 2 ) then

    do j = 1, segment_length
      angle = 2.0D+00 * pi &
        * real ( segment_length - j, kind = 8 ) &
        / real ( segment_length - 1, kind = 8 )
      segment(1,j) = center(1) + r2 * cos ( angle )
      segment(2,j) = center(2) + r2 * sin ( angle )
    end do

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P02_BOUNDARY_SEGMENT - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal SEGMENT_INDEX = ', segment_index
    stop

  end if

  return
end
subroutine p02_boundary_segment_length ( segment_index, h, segment_length )

!*****************************************************************************80
!
!! P02_BOUNDARY_SEGMENT_LENGTH returns boundary segment lengths in problem 02.
!
!  Modified:
!
!    13 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer SEGMENT_INDEX, the index of one of the boundary segments.
!
!    Input, real ( kind = 8 ) H, the suggested spacing between points.
!
!    Output, integer SEGMENT_LENGTH, the number of points in the segment.
!
  implicit none

  real ( kind = 8 ) h
  integer n
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  real ( kind = 8 ), parameter :: r2 = 0.4D+00
  integer segment_index
  integer segment_length

  if ( h <= 0.0D+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P02_BOUNDARY_SEGMENT_LENGTH - Fatal error!'
    write ( *, '(a,g14.6)' ) '  Nonpositive H = ', h
    stop
  end if

  if ( segment_index == 1 ) then

    n = nint ( 2.0D+00 * pi * r1 / h )
    n = max ( n, 5 )
    segment_length = n

  else if ( segment_index == 2 ) then

    n = nint ( 2.0D+00 * pi * r2 / h )
    n = max ( n, 5 )
    segment_length = n
    
  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P02_BOUNDARY_SEGMENT_LENGTH - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal SEGMENT_INDEX = ', segment_index
    stop
  end if

  return
end
subroutine p02_boundary_segment_num ( boundary_segment_num )

!*****************************************************************************80
!
!! P02_BOUNDARY_SEGMENT_NUM counts the boundary segments in problem 02.
!
!  Modified:
!
!    15 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Output, integer BOUNDARY_SEGMENT_NUM, the number of boundary segments.
!
  implicit none

  integer boundary_segment_num

  boundary_segment_num = 2

  return
end
subroutine p02_box ( m, lo, hi )

!*****************************************************************************80
!
!! P02_BOX returns a bounding box for problem 02.
!
!  Modified:
!
!    15 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Output, real ( kind = 8 ) LO(M), HI(M), coordinates of the
!    low and high corners of the box.
!
  implicit none

  integer m

  real ( kind = 8 ), dimension ( 2 ) :: center = (/ 0.0D+00, 0.0D+00 /)
  real ( kind = 8 ) hi(m)
  real ( kind = 8 ) lo(m)
  real ( kind = 8 ), parameter :: r1 = 1.0D+00
  real ( kind = 8 ), parameter :: r2 = 0.4D+00

  lo(1:m) = (/ center(1) - r1, center(2) - r1 /)
  hi(1:m) = (/ center(1) + r1, center(2) + r1 /)

  return
end
subroutine p02_density ( m, n, point, density )

!*****************************************************************************80
!
!! P02_DENSITY returns the density for problem 02.
!
!  Modified:
!
!    15 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = 8 ) POINT(M,N), the coordinates 
!    of the points.
!
!    Output, real ( kind = 8 ) DENSITY(N), the mesh density at 
!    each point.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ) density(n)
  real ( kind = 8 ) point(m,n)

  density(1:n) = 1.0D+00

  return
end
subroutine p02_element_size ( element_size )

!*****************************************************************************80
!
!! P02_ELEMENT_SIZE returns a typical element size for problem 02.
!
!  Modified:
!
!    15 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ELEMENT_SIZE, a typical element size.
!
  implicit none

  real ( kind = 8 ) element_size

  element_size = 0.1D+00

  return
end
subroutine p02_fixed_num ( fixed_num )

!*****************************************************************************80
!
!! P02_FIXED_NUM returns the number of fixed points in problem 02.
!
!  Modified:
!
!    15 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Output, integer FIXED_NUM, the number of fixed points.
!
  implicit none

  integer fixed_num

  fixed_num = 0

  return
end
subroutine p02_fixed_points ( m, fixed_num, fixed )

!*****************************************************************************80
!
!! P02_FIXED_POINTS returns the fixed points in problem 02.
!
!  Modified:
!
!    15 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Per-Olof Persson and Gilbert Strang,
!    A Simple Mesh Generator in MATLAB,
!    SIAM Review,
!    Volume 46, Number 2, June 2004, pages 329-345.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer FIXED_NUM, the number of fixed points.
!
!    Output, real ( kind = 8 ) FIXED(M,FIXED_NUM), the coordinates 
!    of the fixed points.
!
  implicit none

  integer m
  integer fixed_num

  real ( kind = 8 ) fixed(m,fixed_num)

  return
end
subroutine p02_header (