program cvt_dense_prb
!
!******************************************************************************
!
!! CVT_DENSE_PRB tests the CVT_DENSE routines.
!
!
!  Modified:
!
!    08 October 2001
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer, parameter :: dim_num = 1
  integer, parameter :: cell_num = 3
  integer, parameter :: sample_num = 10000
!
  real box_max(dim_num)
  real box_min(dim_num)
  integer cell
  real cell_gen_coord(dim_num,cell_num)
  real cell_volume_desired(cell_num)
  real cell_volume(cell_num)
  integer count(cell_num)
  real density_coord(dim_num,cell_num)
  real density_value(cell_num)
  real density_value_max
  integer i
  integer density_it
  integer, parameter :: density_it_max = 10
  integer nearest
  integer ngen
  real, parameter :: region_volume = 10.0
  real sample_coord(dim_num,sample_num)
  integer seed
!
  call timestamp()

  seed = 1952
  call random_initialize ( seed )

  box_min(1) = 0.0
  box_max(1) = 10.0

  cell_volume_desired(1:cell_num) = &
    (/ 2.0, 3.0, 5.0 /)

  cell_volume_desired(1:cell_num) = cell_volume_desired(1:cell_num) &
    * region_volume / sum ( cell_volume_desired(1:cell_num) )
!
!  Set the initial generators randomly.
!
! call region_sampler ( dim_num, cell_num, box_min, box_max, &
!   cell_gen_coord, ngen )
!
!  SORT THEM.
!
 
  cell_gen_coord(1,1:cell_num) = &
    (/ 1.0, 4.0, 8.0 /)

  density_it = 0
  write ( *, * ) ' '
  write ( *, * ) 'Density iteration number: ', density_it 
  write ( *, * ) ' '
  write ( *, * ) 'Cell, Generator'
  write ( *, * ) ' '
  do cell = 1, cell_num
    write ( *, '(i3,3f10.4)' ) cell, cell_gen_coord(1:dim_num,cell)
  end do
!
!  Determine the cell volumes.
!
  call cell_volume_computation ( dim_num, box_min, box_max, cell_num, &
    cell_gen_coord, sample_num, cell_volume, region_volume )

  write ( *, * ) ' '
  write ( *, * ) 'Cell, CELL VOLUME, CELL_VOLUME_DESIRED'
  write ( *, * ) ' '
  do cell = 1, cell_num
    write ( *, '(i4,3f10.4)' ) &
      cell, cell_volume(cell), cell_volume_desired(cell)
  end do
!
!  Initialize the density function.
!
  call density_set ( cell_num, dim_num, box_min, box_max, cell_gen_coord, &
    cell_volume_desired, density_coord, density_value )

  do density_it = 1, density_it_max

    write ( *, * ) ' '
    write ( *, * ) 'Density iteration number: ', density_it 
!
!  Find the CVT, given the current density function.
!
    call cvt_dense_iteration ( dim_num, box_min, box_max, cell_num, &
      cell_gen_coord, density_value, density_coord, sample_num )

    write ( *, * ) ' '
    write ( *, * ) 'Cell, Cell_Generator'
    write ( *, * ) ' '
    do cell = 1, cell_num
      write ( *, '(i3,3f10.4)' ) cell, cell_gen_coord(1:dim_num,cell)
    end do
!
!  Determine the cell volumes.
!
    call cell_volume_computation ( dim_num, box_min, box_max, cell_num, &
      cell_gen_coord, sample_num, cell_volume, region_volume )

    write ( *, * ) ' '
    write ( *, * ) 'Cell, CELL VOLUME, CELL_VOLUME_DESIRED'
    write ( *, * ) ' '
    do cell = 1, cell_num
      write ( *, '(i4,3f10.4)' ) &
        cell, cell_volume(cell), cell_volume_desired(cell)
    end do
!
!  Update the density function.
!
    call density_set ( cell_num, dim_num, box_min, box_max, cell_gen_coord, &
      cell_volume_desired, density_coord, density_value )

  end do

  stop
end
subroutine density_set ( cell_num, dim_num, box_min, box_max, cell_gen_coord, &
  cell_volume_desired, density_coord, density_value )
!
!******************************************************************************
!
!! DENSITY_SET updates the density function.
!
!
!  Discussion:
!
!    In 1D cases, the volume of a cell tends to be proportional to
!    the integral over the cell of 1/Density**3.
!
!    If a cell has a smaller area than we wish, we can try DECREASING
!    the density function locally.
!
!  Modified:
!
!    03 November 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real box_max(cell_num)
  real box_min(cell_num)
  integer cell
  real cell_gen_coord(dim_num,cell_num)
  real cell_volume_desired(cell_num)
  logical, parameter :: debug = .true.
  real density_coord(1:dim_num,cell_num)
  real density_value(cell_num)
  real density_value_max
  integer i
  integer int_num
  real int_x(cell_num+1)
  real int_v(cell_num)
  real temp
!
!  Compute the interval boundaries.
!
  int_num = cell_num
  int_x(1) = box_min(1)
  do i = 2, int_num
    int_x(i) = 0.5 * ( cell_gen_coord(1,i-1) + cell_gen_coord(1,i) )
  end do
  int_x(cell_num+1) = box_max(1)
!
!  Compute the desired integral values.
!
  int_v(1:cell_num) = ( 1.0 / cell_volume_desired(1:cell_num) )**3

  if ( debug ) then

    write ( *, * ) ' '
    write ( *, * ) 'DENSITY_SET:'
    write ( *, * ) '  Determine the desired integral values.'
    write ( *, * ) ' '
    write ( *, * ) 'Cell, Left_X, Right_X, Desired, Achieved'
    write ( *, * ) ' '

    do cell = 1, cell_num
      temp = ( 1.0 / ( int_x(cell+1) - int_x(cell) ) )**3
      write ( *, '(i3,4f10.4)' ) cell, int_x(cell), int_x(cell+1), &
        int_v(cell), temp
    end do

  end if
!
!  Determine the appropriate linear spline for the density.
!
  call spline_linear_intset ( int_num, int_x, int_v, cell_num, &
    density_coord, density_value )
!
!  Normalize so that the maximum value is 1.
!
  density_value_max = maxval ( density_value(1:cell_num) )
  density_value(1:cell_num) = density_value(1:cell_num) / density_value_max
!
!  Optional printout.
!
  if ( debug ) then
    write ( *, * ) ' '
    write ( *, * ) 'DENSITY_SET:'
    write ( *, * ) '  Compute the corresponding density spline.'
    write ( *, * ) ' '
    write ( *, * ) 'Breakpoint, Density_Coord, Density_Value'
    write ( *, * ) ' '

    do cell = 1, cell_num
      write ( *, '(i3,3f10.4,f10.4)' ) cell, density_coord(1:dim_num,cell), &
        density_value(cell)
    end do

  end if

  return
end