subroutine cvt_park_iteration ( dim_num, box_min, box_max, cell_num, &  
  cell_generator, cell_size )
!
!******************************************************************************
!
!! CVT_PARK_ITERATION takes one step of the CVT parking iteration.
!
!
!  Discussion:
!
!    The routine is given a set of points, called "generators", which
!    define a tessellation of the region into Voronoi cells.  Each point
!    defines a cell.  Each cell, in turn, has a centroid, but it is
!    unlikely that the centroid and the generator coincide.
!
!    Each time this CVT iteration is carried out, an attempt is made
!    to modify the generators in such a way that they are closer and
!    closer to being the centroids of the Voronoi cells they generate.
!
!    Each generator is assigned a particular capacity, and the corresponding
!    number of sample points are generated.  Each sample point is assigned
!    to the nearest available generator, with the understanding that
!    once a generator has reached its capacity, it is no longer available.
!
!    Once the sampling is completed, the location of all the
!    generators is adjusted.  This step should decrease the discrepancy
!    between the generators and the centroids.
!
!  Modified:
!
!    07 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real BOX_MIN(DIM_NUM), BOX_MAX(DIM_NUM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, integer CELL_NUM, the number of Voronoi cells.
!
!    Input/output, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the Voronoi
!    cell generators.  On output, these have been modified.
!
!    Input, integer CELL_SIZE(CELL_NUM), the maximum number of points
!    that can be accepted by each cell generator.  The total of these
!    entries is the number of sample points that will be generated.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real box_max(1:dim_num)
  real box_min(1:dim_num)
  real c1
  real c2
  integer cell_count(cell_num)
  real cell_generator(dim_num,cell_num)
  real cell_generator2(dim_num,cell_num)
  integer cell_size(cell_num)
  real dist
  integer i
  integer j
  integer k
  integer nearest1
  integer nearest2
  integer ngen
  integer sample_num
  real x(dim_num)
!
  cell_generator2(1:dim_num,1:cell_num) = 0.0E+00
  cell_count(1:cell_num) = 0

  sample_num = sum ( cell_size(1:cell_num) )

  do j = 1, sample_num
!
!  Generate a sampling point X.
!
    call region_sampler ( dim_num, box_min, box_max, x, ngen )
!
!  Find the nearest cell generator.
!
    call find_closest_park ( dim_num, x, cell_num, cell_generator, &
      cell_size, cell_count, nearest1, nearest2 )
!
!  Add X to the averaging data for CELL_GENERATOR(*,NEAREST).
!
    if ( nearest2 > 0 ) then

      cell_generator2(1:dim_num,nearest2) = &
        cell_generator2(1:dim_num,nearest2) + x(1:dim_num)

      if ( nearest1 /= nearest2 ) then
        cell_generator2(1:dim_num,nearest1) = &
          cell_generator2(1:dim_num,nearest1) + &
          cell_generator(1:dim_num,nearest1) - x(1:dim_num)
      end if

    else

      write ( *, * ) '?'
      stop

    end if

  end do
!
!  Compute the new generators.
!
  do j = 1, cell_num

    if ( cell_count(j) /= 0 ) then

      cell_generator2(1:dim_num,j) = cell_generator2(1:dim_num,j) &
        / real ( cell_count(j) )

    end if

  end do

    write ( *, * ) ' '
    write ( *, * ) 'Cell generators:'
    write ( *, * ) ' '
    do i = 1, cell_num
      write ( *, '(i4,3f15.6)' ) i, cell_generator(1:dim_num,i)
      write ( *, '(i4,3f15.6)' ) i, cell_generator2(1:dim_num,i)
    end do

  cell_generator(1:dim_num,1:cell_num) = cell_generator2(1:dim_num,1:cell_num)

  return
end
subroutine find_closest ( dim_num, x, cell_num, cell_generator, nearest )
!
!******************************************************************************
!
!! FIND_CLOSEST finds the Voronoi cell generator closest to a point X.
!
!
!  Discussion:
!
!    This routine finds the closest Voronoi cell generator by checking every
!    one.  For problems with many cells, this process can take the bulk
!    of the CPU time.  Other approaches, which group the cell generators into
!    bins, can run faster by a large factor.
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real X(DIM_NUM), the point to be checked.
!
!    Input, integer CELL_NUM, the number of cell generatorrs.
!
!    Input, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the cell generators.
!
!    Output, integer NEAREST, the index of the nearest cell generators.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real cell_generator(dim_num,cell_num)
  real distance
  real dist_sq
  integer i
  integer nearest
  real x(dim_num)
!
  nearest = 0
  distance = huge ( distance )

  do i = 1, cell_num

    dist_sq = sum ( ( cell_generator(1:dim_num,i) - x(1:dim_num) )**2 )

    if ( dist_sq < distance ) then
      distance = dist_sq
      nearest = i
    end if

  end do

  distance = sqrt ( distance )

  return
end
subroutine find_closest_park ( dim_num, x, cell_num, cell_generator, &
  cell_size, cell_count, nearest1, nearest2 )
!
!******************************************************************************
!
!! FIND_CLOSEST_PARK finds the Voronoi cell generator closest to a point X.
!
!
!  Discussion:
!
!    This routine "wants" to find the nearest cell generator to a given
!    point.  However, each cell generator has a maximum capacity, and 
!    once that capacity is reached, the generator will not "accept"
!    any more points.  
!
!  Modified:
!
!    07 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real X(DIM_NUM), the point to be checked.
!
!    Input, integer CELL_NUM, the number of cell generatorrs.
!
!    Input, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the cell generators.
!
!    Input, integer CELL_SIZE(CELL_NUM), the maximum number of points that
!    can be accepted by each cell generator.
!
!    Input/output, integer CELL_COUNT(CELL_NUM), the number of points
!    associated with each cell generator.  On return, the entry
!    for the generator associated with the input point X will have
!    been incremented.
!
!    Output, integer NEAREST, the index of the nearest "available" cell
!    generator.  If no cell generator was available, NEAREST is set
!    to -1.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  integer cell_count(cell_num)
  real cell_generator(dim_num,cell_num)
  integer cell_size(cell_num)
  real distance1
  real distance2
  real dist_sq
  integer i
  integer nearest1
  integer nearest2
  real x(dim_num)
!
  nearest1 = -1
  nearest2 = -1
  distance1 = huge ( distance1 )
  distance2 = huge ( distance2 )

  do i = 1, cell_num

    dist_sq = sum ( ( cell_generator(1:dim_num,i) - x(1:dim_num) )**2 )

    if ( dist_sq < distance1 ) then 
      distance1 = dist_sq
      nearest1 = i
    end if

    if ( dist_sq < distance2 ) then 
      if ( cell_count(i) < cell_size(i) ) then
        distance2 = dist_sq
        nearest2 = i
      end if
    end if

  end do

  if ( nearest2 > 0 ) then
    cell_count(nearest2) = cell_count(nearest2) + 1
  else
    write ( *, * ) '?'
    stop
  end if

  return
end
subroutine find_closest_size ( dim_num, x, cell_num, cell_generator, &
  cell_size, nearest )
!
!******************************************************************************
!
!! FIND_CLOSEST_SIZE finds the Voronoi cell generator closest to a point X.
!
!
!  Discussion:
!
!    This routine finds the closest Voronoi cell generator by checking every
!    one.  For problems with many cells, this process can take the bulk
!    of the CPU time.  Other approaches, which group the cell generators into
!    bins, can run faster by a large factor.
!
!  Modified:
!
!    03 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real X(DIM_NUM), the point to be checked.
!
!    Input, integer CELL_NUM, the number of cell generatorrs.
!
!    Input, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the cell generators.
!
!    Output, integer NEAREST, the index of the nearest cell generators.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real cell_generator(dim_num,cell_num)
  real cell_size(cell_num)
  real distance
  real dist_sq
  integer i
  integer nearest
  real x(dim_num)
!
  nearest = 0
  distance = huge ( distance )

  do i = 1, cell_num

    dist_sq = sum ( ( cell_generator(1:dim_num,i) - x(1:dim_num) )**2 )
!
!  This is sensible only for 2D problems...
!
    dist_sq = sqrt ( dist_sq ) / sqrt ( cell_size(i) )

    if ( dist_sq < distance ) then
      distance = dist_sq
      nearest = i
    end if

  end do

  return
end
subroutine generator_init ( dim_num, box_min, box_max, cell_num, cell_generator )
!
!******************************************************************************
!
!! GENERATOR_INIT initializes the Voronoi cell generators.
!
!
!  Discussion:
!
!    The points initialized here will be used to generate a tessellation
!    of the region into Voronoi cells.  Each generator point defines a
!    cell.  The CVT algorithm will try to modify these initial generators
!    in such a way that they are also the centroids of the cells they generate.
!
!    It is probably better to use Halton points for the centroids than
!    uniform random values, in the sense that the algorithm will probably
!    converge more quickly.
!
!  Modified:
!
!    07 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real BOX_MIN(DIM_NUM), BOX_MAX(DIM_NUM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, integer CELL_NUM, the number of Voronoi cells.
!
!    Output, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the Voronoi cell
!    generators.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real box_max(dim_num)
  real box_min(dim_num)
  real cell_generator(dim_num,cell_num)
  integer i
  integer ngen
!
  do i = 1, cell_num

    call region_sampler ( dim_num, box_min, box_max, cell_generator(1,i), ngen )

  end do

  return
end
subroutine quality ( dim_num, cell_num, cell_moment, cell_volume, region_volume )
!
!*******************************************************************************
!
!! QUALITY computes some quality measures for a set of points in a region.
!
!
!  Discussion:
!
!    The quality measures report on how evenly spread the points are.
!
!  Modified:
!
!    03 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer CELL_NUM, the number of cell generators.
!
!    Input, real CELL_MOMENT(DIM_NUM,DIM_NUM,CELL_NUM), the second moment
!    matrix for each Voronoi cell.
!
!    Input, real CELL_VOLUME(CELL_NUM), the Voronoi cell volumes.
!
!    Input, real REGION_VOLUME, the volume of the region,
!    as input by the user or estimated by VCM.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real cell_det(cell_num)
  real cell_moment(dim_num,dim_num,cell_num)
  real cell_trace(cell_num)
  real cell_volume(cell_num)
  real dd_l1
  real dd_l2
  real dd_linf
  real rmat_det_2d
  real ev
  real ev_l1
  real ev_l2
  real ev_linf
  integer i
  integer j
  real matrix2(2,2)
  real matrix3(3,3)
  real region_volume
  real tr
  real tr_l1
  real tr_l2
  real tr_linf
!
!  Measure 1: the deviation of the cell volumes from the expected cell volume.
!
  ev = region_volume / real ( cell_num )
  ev_linf = maxval ( abs ( ev - cell_volume(1:cell_num) ) )
  ev_l1 = sum ( abs ( ev - cell_volume(1:cell_num) ) )
  ev_l2 = sqrt ( sum ( ( ev - cell_volume(1:cell_num) )**2 ) )

  write ( *, * ) ' '
  write ( *, * ) 'QUALITY'
  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  Measure #1:'
  write ( *, * ) '    ( Cell_Volume - Expected Cell Volume )'
  write ( *, * ) ' '
  write ( *, * ) '  Expected Cell Volume = ', ev
  write ( *, * ) ' '
  write ( *, * ) '  L1 norm =              ', ev_l1
  write ( *, * ) '  L2 norm =              ', ev_l2
  write ( *, * ) '  L1 norm / N =          ', ev_l1 / real ( cell_num )
  write ( *, * ) '  L2 norm / sqrt ( N ) = ', ev_l2 / sqrt ( real ( cell_num ) )
  write ( *, * ) '  L-Inf norm =           ', ev_linf
!
!  Measure 2: the deviation of the traces of the cell second moment matrices
!  from the average.
!
  do i = 1, cell_num
    cell_trace(i) = 0.0E+00
    do j = 1, dim_num
      cell_trace(i) = cell_trace(i) + cell_moment(j,j,i)
    end do
  end do

  tr = sum ( cell_trace(1:cell_num) ) / real ( cell_num )
  tr_linf = maxval ( abs ( tr - cell_trace(1:cell_num) ) )
  tr_l1 = sum ( abs ( tr - cell_trace(1:cell_num) ) )
  tr_l2 = sqrt ( sum ( ( tr - cell_trace(1:cell_num) )**2 ) )

  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  Measure #2:'
  write ( *, * ) '    ( Cell_Trace - Average Cell Trace )'
  write ( *, * ) ' '
  write ( *, * ) '  Average Cell Trace = ', tr
  write ( *, * ) ' '
  write ( *, * ) '  L1 norm =              ', tr_l1
  write ( *, * ) '  L2 norm =              ', tr_l2
  write ( *, * ) '  L1 norm / N =          ', tr_l1 / real ( cell_num )
  write ( *, * ) '  L2 norm / sqrt ( N ) = ', tr_l2 / sqrt ( real ( cell_num ) )
  write ( *, * ) '  L-Inf norm =           ', tr_linf
!
!  Measure 3: the determinant of the deviatoric matrix
!
  do i = 1, cell_num
    matrix2(1:dim_num,1:dim_num) = cell_moment(1:dim_num,1:dim_num,i)
    do j = 1, dim_num
      matrix2(j,j) = matrix2(j,j) - cell_trace(i) / real ( dim_num )
    end do

    cell_det(i) = rmat_det_2d ( matrix2 )

  end do

  dd_linf = maxval ( abs ( cell_det(1:cell_num) ) )
  dd_l1 = sum ( abs ( cell_det(1:cell_num) ) )
  dd_l2 = sqrt ( sum ( ( cell_det(1:cell_num) )**2 ) )

  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  Measure #3:'
  write ( *, * ) '    ( The determinant of the deviatoric matrix )'
  write ( *, * ) ' '
  write ( *, * ) '  L1 norm =              ', dd_l1
  write ( *, * ) '  L2 norm =              ', dd_l2
  write ( *, * ) '  L1 norm / N =          ', dd_l1 / real ( cell_num )
  write ( *, * ) '  L2 norm / sqrt ( N ) = ', dd_l2 / sqrt ( real ( cell_num ) )
  write ( *, * ) '  L-Inf norm =           ', dd_linf

  return
end
subroutine random_initialize ( seed )
!
!*******************************************************************************
!
!! RANDOM_INITIALIZE initializes the FORTRAN 90 random number seed.
!
!
!  Discussion:
!
!    If you don't initialize the random number generator, its behavior
!    is not specified.  If you initialize it simply by:
!
!      call random_seed
!
!    its behavior is not specified.  On the DEC ALPHA, if that's all you
!    do, the same random number sequence is returned.  In order to actually
!    try to scramble up the random number generator a bit, this routine
!    goes through the tedious process of getting the size of the random
!    number seed, making up values based on the current time, and setting
!    the random number seed.
!
!  Modified:
!
!    03 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer SEED.
!    If SEED is zero on input, then you're asking this routine to come up
!    with a seed value, which is returned as output.
!    If SEED is nonzero on input, then you're asking this routine to
!    use the input value of SEED to initialize the random number generator.
!
  integer count
  integer count_max
  integer count_rate
  integer i
  integer seed
  integer, allocatable :: seed_vector(:)
  integer seed_size
  real t
!
!  Initialize the random number seed.
!
  call random_seed ( )
!
!  Determine the size of the random number seed.
!
  call random_seed ( size = seed_size )
!
!  Allocate a seed of the right size.
!
  allocate ( seed_vector(seed_size) )

  if ( seed /= 0 ) then

    write ( *, * ) ' '
    write ( *, * ) 'RANDOM_INITIALIZE'
    write ( *, * ) '  Initialize RANDOM_NUMBER with user SEED = ', seed

  else

    call system_clock ( count, count_rate, count_max )

    seed = count

    write ( *, * ) ' '
    write ( *, * ) 'RANDOM_INITIALIZE'
    write ( *, * ) '  Initialize RANDOM_NUMBER with arbitrary SEED = ', seed

  end if
!
!  Now set the seed.
!
  seed_vector(1:seed_size) = seed

  call random_seed ( put = seed_vector(1:seed_size) )
!
!  Free up the seed space.
!
  deallocate ( seed_vector )
!
!  Call the random number routine a bunch of times.
!
  do i = 1, 100
    call random_number ( harvest = t )
  end do

  return
end
subroutine region_sampler ( dim_num, box_min, box_max, x, ngen )
!
!******************************************************************************
!
!! REGION_SAMPLER returns a sample point in the physical region.
!
!
!  Discussion:
!
!    The calculations are done in DIM_NUM dimensional space.
!
!    The physical region is enclosed in a bounding box.
!
!    A point is chosen in the bounding box by a uniform random
!    number generator.
!
!    If a user-supplied routine determines that this point is
!    within the physical region, this routine returns.  Otherwise,
!    a new random point is chosen.
!
!  Modified:
!
!    07 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real BOX_MIN(DIM_NUM), BOX_MAX(DIM_NUM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Output, real X(DIM_NUM), the random point.
!
!    Output, integer NGEN, the number of points that were generated.
!    This is at least 1, but may be larger if some points were rejected.
!
  integer dim_num
!
  real box_max(dim_num)
  real box_min(dim_num)
  integer ival
  integer ngen
  real r(dim_num)
  real x(dim_num)
!
  ngen = 0

  do

    ngen = ngen + 1

    if ( ngen > 10000 ) then
      write ( *, * ) ' '
      write ( *, * ) 'REGION_SAMPLER - Fatal error!'
      write ( *, * ) '  Generated ', ngen, ' rejected points in a row.'
      write ( *, * ) '  There may be a problem with the geometry definition.'
      write ( *, * ) ' '
      write ( *, * ) '  Current random value is:'
      write ( *, '(3g14.6)' ) r(1:dim_num)
      write ( *, * ) ' '
      write ( *, * ) '  Current sample point is:'
      write ( *, '(3g14.6)' ) x(1:dim_num)
      stop
    end if
!
!  Generate a point at random.
!
    call random_number ( r(1:dim_num) )
!
!  Determine a point in the bounding box.
!
    x(1:dim_num) = ( ( 1.0E+00 - r(1:dim_num) ) * box_min(1:dim_num) &
                            + r(1:dim_num)   * box_max(1:dim_num) )

    call test_region ( x, dim_num, ival )

    if ( ival == 1 ) then
      exit
    end if

  end do

  return
end
function rmat_det_2d ( a )
!
!*******************************************************************************
!
!! RMAT_DET_2D computes the determinant of a 2 by 2 matrix.
!
!
!  Formula:
!
!    The determinant of a 2 by 2 matrix is
!
!      a11 * a22 - a12 * a21.
!
!  Modified:
!
!    01 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real A(2,2), the matrix whose determinant is desired.
!
!    Output, real RMAT_DET_2D, the determinant of the matrix.
!
  real a(2,2)
  real rmat_det_2d
!
  rmat_det_2d = a(1,1) * a(2,2) - a(1,2) * a(2,1)

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

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

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

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

  return
end
subroutine vcm ( dim_num, box_min, box_max, cell_num, cell_generator, cell_size, &
  region_volume, cell_volume, cell_centroid, cell_moment )
!
!*******************************************************************************
!
!! VCM calculates Voronoi cell volumes, centroids and second moments.
!
!
!  Discussion:
!
!    A Monte Carlo sampling is used to estimate the quantities.
!
!  Modified:
!
!    07 September 2001
!
!  Author:
!
!    Janet Peterson
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real BOX_MIN(DIM_NUM), BOX_MAX(DIM_NUM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, integer CELL_NUM, the number of cell generators.
!
!    Input, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the cell generators.
!
!    Input, integer CELL_SIZE(CELL_NUM), the maximum number of points
!    that can be accepted by each cell generator.  The total of these
!    entries is the number of sample points that will be generated.
!
!    Output, real REGION_VOLUME, the volume of the region.
!
!    Output, real CELL_VOLUME(CELL_NUM), the Voronoi cell volumes.
!
!    Output, real CELL_CENTROID(DIM_NUM,CELL_NUM), the Voronoi cell centroids.
!
!    Output, real CELL_MOMENT(DIM_NUM,DIM_NUM,CELL_NUM), the second moment
!    matrix for each Voronoi cell.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real box_max(dim_num)
  real box_min(dim_num)
  real box_volume
  real cell_centroid(dim_num,cell_num)
  real cell_generator(dim_num,cell_num)
  integer cell_count(cell_num)
  real cell_moment(dim_num,dim_num,cell_num)
  integer cell_size(cell_num)
  real cell_volume(cell_num)
  integer i
  integer j
  integer k
  integer nearest
  integer ngen
  integer ntries
  real region_volume
  real region_volume_estimate
  integer sample_num
  real x(dim_num)
!
!  Zero out the arrays.
!
  cell_centroid(1:dim_num,1:cell_num) = 0.0E+00
  cell_count(1:cell_num) = 0
  cell_moment(1:dim_num,1:dim_num,1:cell_num) = 0.0E+00
!
!  Sample the region SAMPLE_NUM times, and keep track of which 
!  cell generator is closest to each sampling point.
!
  ntries = 0

  sample_num = sum ( cell_size(1:cell_num) )

  do k = 1, sample_num

    call region_sampler ( dim_num, box_min, box_max, x, ngen )
    ntries = ntries + ngen

    call find_closest ( dim_num, x, cell_num, cell_generator, nearest )

    cell_count(nearest) = cell_count(nearest) + 1
    cell_centroid(1:dim_num,nearest) = cell_centroid(1:dim_num,nearest) + x(1:dim_num)

    do i = 1, dim_num
      do j = 1, dim_num
        cell_moment(i,j,nearest) = cell_moment(i,j,nearest) + x(i) * x(j)
      end do
    end do

  end do
!
!  Estimate the area or volume if it was not given.
!
  box_volume = product ( ( box_max(1:dim_num)-box_min(1:dim_num) ) )
  region_volume = real ( sample_num ) * box_volume / real ( ntries )

  write ( *, * ) ' '
  write ( *, * ) '  Volume of bounding box is     ', box_volume
  write ( *, * ) '  Estimated volume of region is ', region_volume
!
!  Estimate the geometric integrals for each Voronoi cell.
!
  do k = 1, cell_num

    if ( cell_count(k) > 0 ) then

      cell_volume(k) = real ( cell_count(k) ) * region_volume &
        / real ( sample_num )
      cell_centroid(1:dim_num,k) = cell_centroid(1:dim_num,k) / real ( cell_count(k) )
      cell_moment(1:dim_num,1:dim_num,k) = cell_moment(1:dim_num,1:dim_num,k) &
        / real ( cell_count(k) )

      do i = 1, dim_num
        do j = 1, dim_num
          cell_moment(i,j,k) = cell_moment(i,j,k) &
            - cell_centroid(i,k) * cell_centroid(j,k)
        end do
      end do

    end if

  end do

  return
end