program cvt_size_prb
!
!******************************************************************************
!
!! CVT_SIZE_PRB runs a test of CVT_SIZE, the sized CVT library.
!
  implicit none
!
! integer, parameter :: cell_num = 2
  integer, parameter :: cell_num = 10
  integer, parameter :: dim_num = 2
!
  real a
  integer area_it
  integer, parameter :: area_it_max = 20
  real area_disc
  real angle
  real b
  real blue
  real, parameter, dimension ( dim_num ) :: box_min = &
    (/ 0.0E+00, 0.0E+00 /)
  real, parameter, dimension ( dim_num ) :: box_max = &
    (/ 10.0E+00, 10.0E+00 /)
  real c
  real cell_generator(dim_num,cell_num)
  integer, parameter :: cell_sample = 500
  real cell_volume(cell_num)
  real cell_volume_desired(cell_num)
  real cell_weight(cell_num)
  integer centroid_it
  integer, parameter :: centroid_it_max = 20
  real factor
  character ( len = 80 ) file_name
  real green
  integer i
  integer ierror
  integer ios
  integer it
  integer iunit
  integer j
  integer, parameter :: maxit = 1
  integer nearest
  integer ngen
  logical, parameter :: plot = .true.
  real red
  real region_volume
  integer sample_num
  integer, save :: seed = 0
  logical, parameter :: write_output = .true.
  real xmax
  real xmin
  real xy(2)
  real ymax
  real ymin
!
!  Print introduction and options.
!
  call timestamp ( )

  region_volume = 100.0
!
!  Set the desired cell volumes.
!
  do i = 1, cell_num
    cell_weight(i) = real ( i )
!   if ( i <= 3 * cell_num / 4 ) then
!     cell_weight(i) = 1.0
!   else
!     cell_weight(i) = 4.0
!   end if
  end do

  cell_weight(1:cell_num) = cell_weight(1:cell_num) &
    / sum ( cell_weight(1:cell_num) )
!
!  We assume that the initial cell weights are in fact the desired
!  relative areas, and use this to set the desired absolute areas.
!
  cell_volume_desired(1:cell_num) = region_volume * cell_weight(1:cell_num)
!
!  Now we take the DIM_NUM-th root of CELL_WEIGHT to get a quantity
!  appropriate for treating linear distances.
!
  cell_weight(1:cell_num) = cell_weight(1:cell_num)**( 1.0E+00 / real ( dim_num ) )
  cell_weight(1:cell_num) = cell_weight(1:cell_num) &
    / sum ( cell_weight(1:cell_num) )

  sample_num = cell_num * cell_sample

  write ( *, * ) ' '
  write ( *, * ) 'CVT_SIZE_PRB'
  write ( *, * ) '  A sample problem for the probabilistic'
  write ( *, * ) '  sized Centroidal Voronoi Tesselation algorithm.'
  write ( *, * ) ' '
  write ( *, * ) '  Given a region in 2D, the problem is to determine'
  write ( *, * ) '  GENERATORS, a set of points which define a division'
  write ( *, * ) '  of the region into Voronoid cells, which are also CENTROIDS'
  write ( *, * ) '  of the Voronoi cells, and which have a certain SIZE.'

  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) 'Geometry parameters:'
  write ( *, * ) '-------------------'
  write ( *, * ) ' '
  write ( *, * ) '  The spatial dimension is DIM_NUM = ', dim_num
  write ( *, * ) ' '
  write ( *, * ) '  The minimum corner of the bounding box is:'
  write ( *, '(2f10.4)' ) box_min(1:dim_num)
  write ( *, * ) '  The maximum corner of the bounding box is:'
  write ( *, '(2f10.4)' ) box_max(1:dim_num)

  write ( *, * ) ' '
  write ( *, * ) 'CVT Algorithm parameters:'
  write ( *, * ) '-------------------------'
  write ( *, * ) ' '
  write ( *, * ) '  The number of Voronoi cells to generate: ', cell_num
  write ( *, * ) '  Number of iterations to determine CVT: ', maxit
  write ( *, * ) '  Number of sampling points: ', sample_num

  write ( *, * ) '  Voronoi cell generators are initialized by RANDOM_NUMBER.'
  write ( *, * ) ' '
  write ( *, * ) '  The desired CVT cell volumes and initial weights are:'
  write ( *, * ) ' '
  do i = 1, cell_num
    write ( *, '(i2,2g14.6)' ) i, cell_volume_desired(i), cell_weight(i)
  end do

  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) 'Nearest Neighbor Search parameters:'
  write ( *, * ) '-----------------------------------'
  write ( *, * ) ' '
  write ( *, * ) '  The nearest neighbor search is not speeded up.'
  write ( *, * ) '  The nearest neighbor search is done by exhaustion.'

  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) 'Miscellaneous parameters:'
  write ( *, * ) '------------------------'
  write ( *, * ) ' '
  if ( write_output ) then
    write ( *, * ) '  Generator and moment output files WILL be written.'
  else
    write ( *, * ) '  Generator and moment output files will NOT be written.'
  end if

  write ( *, * ) ' '
!
!  Initialize the random number generator.
!
!  If SEED is zero on input, then the routine will make up a seed.
!  If SEED is nonzero, then a reproducible sequence of random numbers
!  defined by SEED will be chosen.
!
  call random_initialize ( seed )
!
!  Initialize the Voronoi cell generators.
!
  write ( *, * ) ' '
  write ( *, * ) '  Initializing the cell generators.'

  call generator_init ( dim_num, box_min, box_max, cell_num, cell_generator )
!
!  Carry out the CVT iteration, which drives the Voronoi generators 
!  and Voronoi centroids closer and closer.
!
  write ( *, * ) ' '
  write ( *, * ) '  Carry out the CVT iteration.'

  do area_it = 1, area_it_max

    write ( *, * ) ' '
    write ( *, * ) '  Area iteration ', area_it

    do centroid_it = 1, centroid_it_max

      call cvt_size_iteration ( dim_num, box_min, box_max, cell_num, &
        cell_generator, cell_weight, sample_num, cell_volume, region_volume )

    end do

    write ( *, * ) ' '
    write ( *, * ) '        Desired, Actual, Linear'
    write ( *, * ) '  Cell  Volume   Volume  Weight'
    write ( *, * ) ' '
    do i = 1, cell_num
      write ( *, '(i6,3g14.6)' ) i, cell_volume_desired(i), cell_volume(i), cell_weight(i)

    end do
!
!  Determine the discrepancies.
!
    area_disc = 0.0
    do i = 1, cell_num
      area_disc = area_disc + abs ( cell_volume_desired(i) - cell_volume(i) )
    end do
    write ( *, * ) ' '
    write ( *, * ) 'Area discrepancy = ', area_disc

    do i = 1, cell_num
!     cell_weight(i) = cell_weight(i) * &
!       ( cell_volume_desired(i) / cell_volume(i) )**( 1.0E+00 / real ( dim_num ) )

      factor = ( cell_volume_desired(i) / cell_volume(i) )
      factor = factor**( 1.0E+00 / real ( dim_num ) )
      factor = min ( factor, 1.25 )
      cell_weight(i) = cell_weight(i) * factor
        
    end do

    cell_weight(1:cell_num) = cell_weight(1:cell_num) &
      / sum ( cell_weight(1:cell_num) )

  end do

!  Print some stuff.
!
  if ( .true. ) then
    write ( *, * ) ' '
    write ( *, * ) 'Cell generators:'
    write ( *, * ) ' '
    do i = 1, cell_num
      write ( *, '(i4,3d15.6)' ) i, cell_generator(1:dim_num,i)
    end do
    write ( *, * ) ' '
    write ( *, * ) 'Cell volumes:'
    write ( *, * ) ' '
    do i = 1, cell_num
      write ( *, '(i4,3d15.6)' ) i, cell_volume(i)
    end do
  end if
!
!  Write generators and moments to files.
!
  if ( write_output ) then

    open ( unit = 1, file = 'cvt_generators.dat', status = 'replace', &
      iostat = ios )
    do i = 1, cell_num
      write ( 1, '(3d15.6)' ) cell_generator(1:dim_num,i)
    end do
    close ( unit = 1 )

    open ( unit = 1, file = 'cvt_volume.dat', status = 'replace', iostat = ios )
    write ( 1, '(d15.6)' ) cell_volume(1:cell_num)
    close ( unit = 1 )

  end if
!
!  Make a Voronoi plot
!
  if ( plot ) then

    iunit = 1
    file_name = 'cvt_size_prb.ps'

    call ps_file_open ( file_name, iunit, ierror )

    if ( ierror /= 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'CVT_SIZE_PRB'
      write ( *, * ) '  File creation error ', ierror
      stop
    end if

    call ps_file_head ( file_name )
!
!  Define the size of the page.
!
    xmin = box_min(1)
    ymin = box_min(2)

    xmax = box_max(1)
    ymax = box_max(2)

    call ps_page_head ( xmin, ymin, xmax, ymax )
!
!  Specify the line color and draw lines.
!
    do i = 1, cell_num

      call random_number ( harvest = angle )
      angle = 360.0 * real ( i - 1 ) / real ( cell_num )
      call hexcol ( angle, red, green, blue )

!     call random_number ( harvest = red )
!     call random_number ( harvest = green )
!     call random_number ( harvest = blue )

      call ps_color_line_set ( red, green, blue )

      do j = 1, cell_sample

        do 

          call region_sampler ( dim_num, box_min, box_max, xy, ngen )

          call find_closest_size ( dim_num, xy, cell_num, cell_generator, &
            cell_weight, nearest )

          if ( nearest == i ) then
            exit
          end if

        end do

        call ps_mark_circle ( xy(1), xy(2) )

      end do

    end do
    red = 0.0E+00
    green = 0.0E+00
    blue = 0.0E+00

    call ps_color_line_set ( red, green, blue )

    do i = 1, cell_num
      call ps_mark_disk ( cell_generator(1,i), cell_generator(2,i) )
    end do

!
!  Close up the page and the file.
!
  call ps_page_tail

  call ps_file_tail

  call ps_file_close ( iunit )

  end if

  write ( *, * ) ' '
  write ( *, * ) 'CVT_SIZE_PRB'
  write ( *, * ) '  Normal end of execution.'

  stop
end
subroutine test_region ( x, dim_num, ival )
!
!*******************************************************************************
!
!! TEST_REGION determines if a point is within the physical region.
!
!
!  Discussion:
!
!    Using a simple routine like this is only appropriate for a simple 
!    region that can be easily defined by user formulas.  
!
!    Computation of the "on-the-boundary" case is not considered important.
!    Only "inside" or "outside" is essential.
!
!  Modified:
!
!    16 April 2001
!
!  Parameters:
!
!    Input, real X(DIM_NUM), the point to be checked.
!
!    Input, integer DIM_NUM, the dimension of the space.
!
!    Output, integer IVAL, indicates the status of the point:
!    -1: the point is on the boundary of the region.
!     0: the point is outside the region.
!    +1: the point is inside the region.
!
  implicit none
!
  integer dim_num
!
  integer ival
  real x(dim_num)
!
  if ( 0.0 <= x(1) .and. x(1) <= 10.0 .and. &
       0.0 <= x(2) .and. x(2) <= 10.0 ) then
    ival = 1
  else
    ival = 0
  end if

  return
end
subroutine hexcol ( angle, r, g, b )
!
!*******************************************************************************
!
!! HEXCOL returns a color on the perimeter of the color hexagon.
!
!
!  Modified:
!
!    12 July 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ANGLE, the angle in the color hexagon.  The sextants are
!    defined by the following points:
!        0 degrees, 1, 0, 0, red;
!       60 degrees, 1, 1, 0, yellow;
!      120 degrees, 0, 1, 0, green;
!      180 degrees, 0, 1, 1, cyan;
!      240 degrees, 0, 0, 1, blue;
!      300 degrees, 1, 0, 1, magenta.
!
!    Output, real R, G, B, RGB specifications for the color that lies
!    at the given angle, on the perimeter of the color hexagon.  One
!    value will be 1, and one value will be 0.
!
  real angle
  real b
  real g
  real r
!
  angle = mod ( angle, 360.0E+00 )

  if ( angle < 0.0E+00 ) then
    angle = angle + 360.0E+00
  end if

  if ( angle <= 60.0E+00 ) then
    r = 1.0E+00
    g = angle / 60.0E+00
    b = 0.0E+00
  else if ( angle <= 120.0E+00 ) then
    r = ( 120.0E+00 - angle ) / 60.0E+00
    g = 1.0E+00
    b = 0.0E+00
  else if ( angle <= 180.0E+00 ) then
    r = 0.0E+00
    g = 1.0E+00
    b = ( angle - 120.0E+00 ) / 60.0E+00
  else if ( angle <= 240.0E+00 ) then
    r = 0.0E+00
    g = ( 240.0E+00 - angle ) / 60.0E+00
    b = 1.0E+00
  else if ( angle <= 300.0E+00 ) then
    r = ( angle - 240.0E+00 ) / 60.0E+00
    g = 0.0E+00
    b = 1.0E+00
  else if ( angle <= 360.0E+00 ) then
    r = 1.0E+00
    g = 0.0E+00
    b = ( 360.0E+00 - angle ) / 60.0E+00
  end if

  return
end