program cvt_park_prb
!
!******************************************************************************
!
!! CVT_PARK_PRB runs a test of CVT_PARK, the parking lot CVT library.
!
  implicit none
!
  integer, parameter :: cell_num = 2
  integer, parameter :: dim_num = 1
!
  real a
  real angle
  real b
  real blue
  real, parameter, dimension ( 2 ) :: box_min = &
    (/ 0.0E+00, 0.0E+00 /)
  real, parameter, dimension ( 2 ) :: box_max = &
    (/ 10.0E+00, 10.0E+00 /)
  real c
  real c_sum
  real cell_centroid(dim_num,cell_num)
  real cell_generator(dim_num,cell_num)
  real cell_moment(cell_num)
  integer cell_size(cell_num)
  integer cell_size_sum
  real cell_volume(cell_num)
  integer clock0
  integer clock1
  integer clock_max
  integer clock_rate
  real etime0
  real etime1
  character ( len = 80 ) file_name
  real green
  integer i
  integer ierror
  integer ios
  integer it
  integer iunit
  integer j
  integer, parameter :: maxit = 4
  integer nearest
  integer ngen
  integer ns
  logical, parameter :: plot = .false.
  logical, save :: quality_checks = .true.
  real red
  real region_volume
  integer, save :: seed = 0
  real tarray(2)
  real time0
  real time1
  logical, parameter :: write_output = .false.
  real xmax
  real xmin
  real xy(2)
  real ymax
  real ymin
!
  call timestamp ( )

  region_volume = 100.0

  ns = 10

  do i = 1, cell_num
    if ( i <= 3 * cell_num / 4 ) then
      cell_size(i) = 1 * ns
    else
      cell_size(i) = 10 * ns
    end if
  end do

  cell_size(1) = 2 * ns
  cell_size(2) = 8 * ns

  cell_size_sum = sum ( cell_size(1:cell_num) )

  write ( *, * ) ' '
  write ( *, * ) 'CVT_PARK_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 ( *, * ) '  Total number of sampling points: ', cell_size_sum

  write ( *, * ) ' '
  write ( *, * ) '  The CVT cell sizes are:'
  write ( *, * ) ' '
  do i = 1, cell_num
    write ( *, * ) i, cell_size(i)
  end do
  write ( *, * ) ' '
  write ( *, * ) '  Cell size sum = ', cell_size_sum
  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 )
!
!  Begin timing.
!
!  The F90 CPU_TIME routine, at least on the DEC_ALPHA, can "wrap around"
!  after a relatively short time, giving negative timings.  So we also
!  call ETIME, which should be measuring the same thing, for checking.
!
!  The ETIME routine is called with an appended underscore, because
!  A) We're compiling with NOUNDERSCORE in order to interface with C;
!  B) but the UFOR library containing ETIME expects a reference to ETIME_.
!
!  The F90 SYSTEM_CLOCK routine should measure real execution time.
!
  call system_clock ( clock0, clock_rate, clock_max )
  call cpu_time ( time0 )
  call etime ( tarray )
  etime0 = tarray(1) + tarray(2)
!
!  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 it = 1, maxit

    call cvt_park_iteration ( dim_num, box_min, box_max, cell_num, &
      cell_generator, cell_size )

    write ( *, * ) ' '
    write ( *, * ) 'Cell generators:'
    write ( *, * ) ' '
    do i = 1, cell_num
      write ( *, '(i4,3f15.6)' ) i, cell_generator(1:dim_num,i)
    end do
!
!  Calculate cell moments.
!
    call vcm ( dim_num, box_min, box_max, cell_num, cell_generator, cell_size, &
      region_volume, cell_volume, cell_centroid, cell_moment )

    write ( *, * ) ' '
    write ( *, * ) '  Desired, Actual, Discrepancy (D-A)/Region:'
    write ( *, * ) ' '
    c_sum = 0.0
    do i = 1, cell_num
      a = real ( cell_size(i) ) * region_volume / real ( cell_size_sum )
      b = cell_volume(i)
      c = abs ( a - b ) / region_volume
      c_sum = c_sum + c
      write ( *, * ) i, a, b, c
    end do
    write ( *, * ) ' '
    write ( *, * ) 'Total discrepancy is ', c_sum

  end do

  stop
!
!  End timing the "interesting" part.
!
! call etime ( tarray )
! etime1 = tarray(1) + tarray(2)
! call cpu_time ( time1 )
! call system_clock ( clock1, clock_rate, clock_max )

! write ( *, * ) ' '
! write ( *, * ) 'Elapsed CPU time, CPU_TIME: ', time1-time0, ' seconds.'
! write ( *, * ) 'Elapsed CPU time, ETIME:    ', etime1-etime0, ' seconds.'
! write ( *, * ) 'Elapsed time, SYSTEM_CLOCK: ', real ( clock1 - clock0 ) &
!   / real ( clock_rate ), ' seconds.'
!
!  Compute quality checks.
!
! if ( quality_checks ) then
!   call quality ( dim_num, cell_num, cell_moment, cell_volume, region_volume )
! end if
!
!  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 )

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

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

    iunit = 1
    file_name = 'cvt_park_prb.ps'

    call ps_file_open ( file_name, iunit, ierror )

    if ( ierror /= 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'CVT_PARK_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.
!
    red = 1.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_circle ( cell_centroid(1,i), cell_centroid(2,i) )
    end do

    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, 2000

        do 

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

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

          if ( nearest == i ) then
            exit
          end if

        end do

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

      end do

    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_PARK_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)
!
  ival = 0

  if ( dim_num == 1 ) then

    if ( 0.0 <= x(1) .and. x(1) <= 10.0 ) then
      ival = 1
    end if

  else if ( dim_num == 2 ) then

    if ( 0.0 <= x(1) .and. x(1) <= 10.0 .and. &
         0.0 <= x(2) .and. x(2) <= 10.0 ) then
      ival = 1
    end if

  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