program mless_main
!
!*******************************************************************************
!
!! MLESS_MAIN is the main program for the meshless basis function routines.
!
  implicit none
!
  double precision, allocatable :: a(:)
  double precision, save, dimension ( 2 ) :: alpha = (/ 0.0D+00, 1.0D+00 /)
  double precision, allocatable :: b(:)
  double precision, allocatable :: basis_center(:,:)
  character ( len = 80 ) :: basis_file = 'basis.dat'
  integer, save :: basis_m = 0
  integer, save :: basis_n = 0
  double precision, allocatable :: basis_radius(:)
  double precision, save :: basis_radius_factor = 1.0D+00
  double precision, save, dimension ( 2 ) :: beta = (/ 0.0D+00, 1.0D+00 /)
  integer, save :: center_choice = 0
  character ( len = 80 ) command
  double precision, allocatable :: cvt_center(:,:)
  character ( len = 80 ) :: cvt_file = 'cvt.dat'
  integer cvt_m
  integer, save :: cvt_maxit = 0
  integer, save :: cvt_n = 0
  integer, save :: cvt_ns = 0
  double precision, allocatable :: cvt_radius(:)
  double precision, save :: cvt_radius_factor = 1.0D+00
  character ( len = 8 ) date
  integer, save :: density_function = 0
  character ( len = 80 ) :: halton_file = 'halton.dat'
  integer, save :: halton_n = 0
  double precision, allocatable :: halton_points(:,:)
  integer i_temp
  integer ierror
  integer ios
  integer last
  integer, save :: maxit = 0
  integer, save :: myrank = 0
  integer, save :: n = 0
  integer ncolumn
  integer, save :: ndim = 2
  integer, save :: ns = 0
  integer, save :: radius_choice = 0
  double precision, save :: radius_factor = 1.0
  double precision, save :: rxmax = 1.0
  double precision, save :: rxmin = -1.0
  double precision, save :: rymax = 1.0
  double precision, save :: rymin = -1.0
  logical s_eqi
  character ( len = 10 ) time
  character ( len = 80 ) :: uniform_file = 'uniform.dat'
  integer, save :: uniform_n = 0
  double precision, allocatable :: uniform_points(:,:)
  character ( len = 80 ) what
!
  call date_and_time ( date, time )

  write ( *, * ) ' '
  write ( *, * ) 'MLESS'
  write ( *, * ) '  Centroidal Voronoi Computations'
  write ( *, * ) ' '
  write ( *, * ) '  Today''s date: ', date
  write ( *, * ) '  Today''s time: ', time
!
!  Initialize the random number generator.
!
  call set_random_seed ( myrank )

  do

    write ( *, * ) ' '
    write ( *, * ) 'Enter command (or HELP for a list)'

    read ( *, '(a)', iostat = ios ) command

    if ( ios /= 0 ) then
      exit
    end if

    if ( s_eqi ( command, 'basis_make' ) ) then

      write ( *, * ) ' '
      write ( *, * ) 'Enter the number of basis functions to create:'

      read ( *, *, iostat = ios ) basis_n

      if ( ios /= 0 ) then
        exit
      end if

      if ( allocated ( basis_center ) ) then
        deallocate ( basis_center )
      end if

      if ( allocated ( basis_radius ) ) then
        deallocate ( basis_radius )
      end if

      allocate ( basis_center(ndim,basis_n) )
      allocate ( basis_radius(basis_n) )

      basis_m = 7 * int ( sqrt ( dble ( basis_n ) ) )

      write ( *, * ) ' '
      write ( *, * ) 'Mysterious parameter BASIS_M = ', basis_m

      do

        write ( *, * ) ' '
        write ( *, * ) 'Initial center point determination:'
        write ( *, * ) ' '
        write ( *, * ) '  1: use a uniform distribution;'
        write ( *, * ) '  2: use a Halton distribution.'
        read ( *, * ) center_choice

        if ( center_choice == 1 ) then
          call uniform_make ( ndim, a, b, basis_n, density_function, &
            basis_center )
          exit
        else if ( center_choice == 2 ) then
          call halton_make ( ndim, basis_n, basis_center )
          exit
        end if

      end do

      do

        write ( *, * ) ' '
        write ( *, * ) 'Initial radius determination:'
        write ( *, * ) ' '
        write ( *, * ) '  1: radius method 1;'
        write ( *, * ) '  2: radius method 2.'

        read ( *, * ) radius_choice

        if ( radius_choice == 1 ) then
          call radius_make_1 ( ndim, a, b, basis_n, basis_m, density_function, &
            basis_center, basis_radius )
          exit
        else if ( radius_choice == 2 ) then
          call radius_make_2 ( ndim, a, b, basis_n, basis_m, basis_center, &
            basis_radius )
          exit
        end if

      end do

    else if ( s_eqi ( command, 'basis_overlap' ) ) then

      call basis_overlap ( ndim, basis_n, basis_center, basis_radius )

    else if ( s_eqi ( command, 'basis_plot' ) ) then

      call basis_plot ( ndim, basis_n, basis_center, basis_radius, a, b )

    else if ( s_eqi ( command, 'basis_read' ) ) then

      call file_column_count ( basis_file, ncolumn )

      ndim = ncolumn - 1

      write ( *, * ) ' '
      write ( *, * ) 'Setting NDIM to ', ndim

      call file_line_count ( basis_file, basis_n )

      write ( *, * ) ' '
      write ( *, * ) 'Setting number of center points to ', basis_n

      if ( allocated ( basis_center ) ) then
        deallocate ( basis_center ) 
      end if
      allocate ( basis_center ( ndim, basis_n ) )

      if ( allocated ( basis_radius ) ) then
        deallocate ( basis_radius ) 
      end if
      allocate ( basis_radius ( basis_n ) )

      call basis_read ( ndim, basis_n, basis_center, basis_radius, basis_file )

      if ( allocated ( a ) ) then
        deallocate ( a ) 
      end if
      allocate ( a ( ndim ) )

      if ( allocated ( b ) ) then
        deallocate ( b ) 
      end if
      allocate ( b ( ndim ) )

      call basis_box ( ndim, basis_n, basis_center, basis_radius, a, b )

    else if ( s_eqi ( command, 'basis_scale' ) ) then

      write ( *, * ) ' '
      write ( *, * ) 'Enter a scale factor to apply to the radius'
      write ( *, * ) 'of the current basis function set:'
      read ( *, *, iostat = ios ) basis_radius_factor
      if ( ios /= 0 ) then
        exit
      end if

      basis_radius(1:basis_n) = basis_radius_factor * basis_radius(1:basis_n)

    else if ( s_eqi ( command, 'basis_write' ) ) then

      if ( basis_n <= 0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'Use the command radius_make first!'
        cycle
      end if

      call basis_write ( ndim, basis_n, basis_center, basis_radius, basis_file )

    else if ( s_eqi ( command, 'center_plot' ) ) then

      call center_plot ( ndim, basis_n, basis_center, basis_radius, a, b )

    else if ( s_eqi ( command, 'cvt_make' ) ) then

      write ( *, * ) ' '
      write ( *, * ) 'Enter the number of CVT points to create:'

      read ( *, *, iostat = ios ) cvt_n
      if ( ios /= 0 ) then
        exit
      end if

      if ( allocated ( a ) ) then
        deallocate ( a )
      end if
      allocate ( a(ndim) )
      a(1:ndim) = -1.0D+00

      if ( allocated ( b ) ) then
        deallocate ( b )
      end if
      allocate ( b(ndim) )
      b(1:ndim) = +1.0D+00

      if ( allocated ( cvt_center ) ) then
        deallocate ( cvt_center )
      end if

      allocate ( cvt_center(ndim,cvt_n) )

      if ( allocated ( cvt_radius ) ) then
        deallocate ( cvt_radius )
      end if

      allocate ( cvt_radius(cvt_n) )

      cvt_m = 7 * int ( sqrt ( dble ( cvt_n ) ) )

      write ( *, * ) ' '
      write ( *, * ) 'Mysterious parameter CVT_M = ', cvt_m

      do

        write ( *, * ) ' '
        write ( *, * ) 'Initial center point determination:'
        write ( *, * ) ' '
        write ( *, * ) '  1: use a uniform distribution;'
        write ( *, * ) '  2: use a Halton distribution.'
        read ( *, * ) center_choice

        if ( center_choice == 1 ) then
          call uniform_make ( ndim, a, b, cvt_n, density_function, cvt_center )
          exit
        else if ( center_choice == 2 ) then
          call halton_make ( ndim, cvt_n, cvt_center )
          exit
        end if

      end do

      do

        write ( *, * ) ' '
        write ( *, * ) 'Initial radius determination:'
        write ( *, * ) ' '
        write ( *, * ) '  1: radius method 1;'
        write ( *, * ) '  2: radius method 2.'

        read ( *, * ) radius_choice

        if ( radius_choice == 1 ) then
          call radius_make_1 ( ndim, a, b, cvt_n, cvt_m, density_function, &
            cvt_center, cvt_radius )
          exit
        else if ( radius_choice == 2 ) then
          call radius_make_2 ( ndim, a, b, cvt_n, cvt_m, cvt_center, cvt_radius )
          exit
        end if

      end do

      call cvt_make ( alpha, beta, ndim, a, b, maxit, ns, density_function, &
        cvt_n, cvt_center )

    else if ( s_eqi ( command, 'cvt_overlap' ) ) then

      call basis_overlap ( ndim, cvt_n, cvt_center, cvt_radius )

    else if ( s_eqi ( command, 'cvt_plot' ) ) then

      call basis_plot ( ndim, cvt_n, cvt_center, cvt_radius, a, b )

    else if ( s_eqi ( command, 'cvt_read' ) ) then

      call cvt_read ( ndim, cvt_n, cvt_center, cvt_radius, cvt_file )

    else if ( s_eqi ( command, 'cvt_scale' ) ) then

      write ( *, * ) ' '
      write ( *, * ) 'Enter a scale factor to apply to the radius'
      write ( *, * ) 'of the current CVT function set:'
      read ( *, *, iostat = ios ) cvt_radius_factor
      if ( ios /= 0 ) then
        exit
      end if

      cvt_radius(1:cvt_n) = cvt_radius_factor * cvt_radius(1:cvt_n)

    else if ( s_eqi ( command, 'cvt_write' ) ) then

      if ( cvt_n <= 0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'Use the command cvt_make first!'
        cycle
      end if

      call cvt_write ( ndim, cvt_n, cvt_center, cvt_radius, cvt_file )

    else if ( s_eqi ( command(1:17), 'DENSITY_FUNCTION=' ) ) then

      call s_to_i ( command(18:), i_temp, ierror, last )

      if ( ierror == 0 ) then

        density_function = i_temp

        write ( *, * ) ' '
        write ( *, * ) '  DENSITY_FUNCTION has been set to ', &
          density_function

      end if

    else if ( s_eqi ( command, 'halton_make' ) ) then

      write ( *, * ) ' '
      write ( *, * ) 'Enter the number of Halton points to create:'
      read ( *, *, iostat = ios ) halton_n
      if ( ios /= 0 ) then
        exit
      end if

      if ( allocated ( halton_points ) ) then
        deallocate ( halton_points )
      end if

      allocate ( halton_points(ndim,halton_n) )

      call halton_make ( ndim, halton_n, halton_points )

    else if ( s_eqi ( command, 'halton_read' ) ) then

      call halton_read ( ndim, halton_n, halton_points, halton_file )

    else if ( s_eqi ( command, 'halton_write' ) ) then

      if ( halton_n <= 0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'Use the command halton_make first!'
        cycle
      end if

      call halton_write ( ndim, halton_n, halton_points, halton_file )

    else if ( s_eqi ( command, 'HELP' ) ) then

      write ( *, * ) ' '
      write ( *, * ) 'Commands:'
      write ( *, * ) ' '
      write ( *, * ) '  DENSITY_FUNCTION = ... (specify density)'
      write ( *, * ) '  MAXIT = ... (specify number of iterations)'
      write ( *, * ) '  N =     ... (specify number of points)'
      write ( *, * ) '  NDIM =  ... (specify spatial dimension)'
      write ( *, * ) '  NS =    ... (specify I don''t know what)'
      write ( *, * ) ' '
      write ( *, * ) '  basis_make'
      write ( *, * ) '  basis_overlap (Count basis function overlaps)'
      write ( *, * ) '  basis_plot'
      write ( *, * ) '  basis_read'
      write ( *, * ) '  basis_scale (make basis functions "wider")'
      write ( *, * ) '  basis_write'
      write ( *, * ) ' '
      write ( *, * ) '  center_plot'
      write ( *, * ) ' '
      write ( *, * ) '  cvt_make'
      write ( *, * ) '  cvt_overlap (Count basis function overlaps)'
      write ( *, * ) '  cvt_plot'
      write ( *, * ) '  cvt_read'
      write ( *, * ) '  cvt_scale (make cvt functions "wider")'
      write ( *, * ) '  cvt_write'
      write ( *, * ) ' '
      write ( *, * ) '  halton_make'
      write ( *, * ) '  halton_read'
      write ( *, * ) '  halton_write'
      write ( *, * ) ' '
      write ( *, * ) '  uniform_make'
      write ( *, * ) '  uniform_read'
      write ( *, * ) '  uniform_write'
      write ( *, * ) ' '
      write ( *, * ) '  Help'
      write ( *, * ) '  Print'
      write ( *, * ) '  Quit'

    else if ( s_eqi ( command(1:6), 'MAXIT=' ) ) then

      call s_to_i ( command(7:), i_temp, ierror, last )

      if ( ierror == 0 ) then

        maxit = i_temp
        write ( *, * ) '  MAXIT has been set to ', maxit

      end if

    else if ( s_eqi ( command(1:2), 'N=' ) ) then

      call s_to_i ( command(3:), i_temp, ierror, last )

      if ( ierror == 0 ) then

        n = i_temp
        write ( *, * ) '  N has been set to ', n

      end if

    else if ( s_eqi ( command(1:5), 'NDIM=' ) ) then

      call s_to_i ( command(6:), i_temp, ierror, last )

      if ( ierror == 0 ) then
        ndim = i_temp
        write ( *, * ) '  NDIM has been set to ', ndim
      end if

    else if ( s_eqi ( command(1:3), 'NS=' ) ) then

      call s_to_i ( command(4:), i_temp, ierror, last )

      if ( ierror == 0 ) then

        ns = i_temp
        write ( *, * ) '  NS has been set to ', ns

      end if

    else if ( s_eqi ( command(1:5), 'Print' ) ) then

      what = adjustl ( command(6:) )
      if ( len_trim ( what ) == 0 ) then
        what = 'ALL'
      end if

      if ( s_eqi ( what, 'N' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, * ) '  The number of points, N = ', n
      end if

      if ( s_eqi ( what, 'NDIM' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, * ) '  The spatial dimension, NDIM = ', ndim
      end if

      if ( s_eqi ( what, 'MAXIT' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, * ) '  The maximum number of iterations, MAXIT = ', maxit
      end if

      if ( s_eqi ( what, 'NS' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, * ) '  I don''t know what this is, NS = ', ns
      end if

      if ( s_eqi ( what, 'DENSITY_FUNCTION' ) .or. s_eqi ( what, 'ALL' ) ) then
        write ( *, * ) '  The density function, DENSITY_FUNCTION = ', &
          density_function
      end if

    else if ( s_eqi ( command, 'Quit' ) ) then

      exit

    else if ( s_eqi ( command, 'uniform_make' ) ) then

      write ( *, * ) ' '
      write ( *, * ) 'Enter the number of uniform points to create:'

      read ( *, *, iostat = ios ) uniform_n

      if ( ios /= 0 ) then
        exit
      end if

      if ( allocated ( uniform_points ) ) then
        deallocate ( uniform_points )
      end if

      allocate ( uniform_points(ndim,uniform_n) )

      call uniform_make ( ndim, a, b, uniform_n, density_function, &
        uniform_points )

    else if ( s_eqi ( command, 'uniform_read' ) ) then

      call uniform_read ( ndim, uniform_n, uniform_points, uniform_file )

    else if ( s_eqi ( command, 'uniform_write' ) ) then

      if ( uniform_n <= 0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'Use the command uniform_make first!'
        cycle
      end if

      call uniform_write ( ndim, uniform_n, uniform_points, uniform_file )

    else

      write ( *, * ) ' '
      write ( *, * ) 'MLESS - Warning'
      write ( *, * ) '  Your command was not recognized.'

    end if

  end do

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

  stop
end
subroutine basis_box ( ndim, n, center, radius, a, b )
!
!*******************************************************************************
!
!! BASIS_BOX finds a box that contains all the basis functions.
!
!
!  Modified:
!
!    26 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, double precision CENTER(NDIM,N), the basis function centers.
!
!    Input, double precision RADIUS(N), the basis function radii.
!
!    Output, double precision A(NDIM), B(NDIM), the minimum and maximum
!    values in each dimension.
!
  implicit none
!
  integer n
  integer ndim
!
  double precision a(ndim)
  double precision b(ndim)
  double precision, dimension ( ndim, n ) :: center
  integer i
  double precision radius(n)
!
  do i = 1, ndim
    a(i) = minval ( center(i,1:n) - radius(1:n) )
    b(i) = maxval ( center(i,1:n) + radius(1:n) )
  end do

  write ( *, * ) ' '
  write ( *, * ) 'BASIS_BOX'
  write ( *, * ) '  Compute bounding box.'
  write ( *, * ) ' '
  write ( *, * ) '  Box limits for center points:'
  write ( *, * ) ' '
  write ( *, * ) '  Dimension  Lower  Upper'
  write ( *, * ) ' '
  do i = 1, ndim
    write ( *, * ) i, minval ( center(i,1:n) ), maxval ( center(i,1:n) )
  end do

  write ( *, * ) ' '
  write ( *, * ) '  Box limits for radius:'
  write ( *, * ) ' '
  write ( *, * ) '  Lower  Upper'
  write ( *, * ) ' '
  write ( *, * ) minval ( radius(1:n) ), maxval ( radius(1:n) )

  write ( *, * ) ' '
  write ( *, * ) '  Box limits for basis support:'
  write ( *, * ) ' '
  write ( *, * ) '  Dimension  Lower  Upper'
  write ( *, * ) ' '

  do i = 1, ndim
    write ( *, * ) i, a(i), b(i)
  end do

  return
end
subroutine basis_overlap ( ndim, n, center, radius )
!
!*******************************************************************************
!
!! BASIS_OVERLAP counts the number of basis functions with overlapping support.
!
!
!  Modified:
!
!    19 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of basis functions.
!
!    Input, double precision CENTER(NDIM,N), the basis function centers.
!
!    Input, double precision RADIUS(N), the basis function radii.
!
  implicit none
!
  integer n
  integer ndim
!
  double precision, dimension ( ndim, n ) :: center
  double precision dist
  integer i
  integer j
  integer mel
  integer nel
  double precision, dimension ( n ) :: radius
!
  if ( n <= 1 ) then
    write ( *, * ) ' '
    write ( *, * ) 'BASIS_OVERLAP - Warning!'
    write ( *, * ) '  N <= 1, so no overlap possible.'
    return
  end if

  mel = ( n * ( n - 1 ) ) / 2
  nel = 0

  do i = 1, n
    do j = i+1, n

      dist = sqrt ( sum ( ( center(1:ndim,i) - center(1:ndim,j) )**2 ) )

      if ( dist <= ( radius(i) + radius(j) ) ) then
        nel = nel + 1
      end if

    end do
  end do

  write ( *, * ) ' '
  write ( *, * ) 'BASIS_OVERLAP:'
  write ( *, * ) '  Count the pairs of basis functions with overlapping support.'
  write ( *, * ) ' '
  write ( *, * ) '  Maximum possible count = ', mel
  write ( *, * ) '  Actual count =           ', nel
  write ( *, * ) '  Percentage =             ', real ( 100 * nel ) / real ( mel )

  return
end
subroutine basis_plot ( ndim, n, center, radius, a, b )
!
!*******************************************************************************
!
!! BASIS_PLOT plots the basis functions in the region.
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of basis functions.
!
!    Input, double precision CENTER(NDIM,N), the basis function centers.
!
!    Input, double precision RADIUS(N), the basis function radii.
!
!    Input, double precision A(NDIM), B(NDIM), the minimum and maximum
!    values in each dimension.
!
  implicit none
!
  integer n
  integer ndim
!
  double precision a(ndim)
  double precision b(ndim)
  real blue
  double precision center(ndim,n)
  character ( len = 80 ) file_name
  logical, parameter :: filled = .false.
  real green
  integer i
  integer ierror
  integer iunit
  integer nx
  integer ny
  real r
  double precision radius(n)
  real red
  double precision rxmax
  double precision rxmin
  double precision rymax
  double precision rymin
  real x
  real xmax
  real xmin
  real y
  real ymax
  real ymin
!
!  NOT CORRECT FOR NDIM > 2
!
  rxmin = a(1)
  rxmax = b(1)
  rymin = a(2)
  rymax = b(2)

  iunit = 1
  file_name = 'basis_plot.ps'

  call ps_file_open ( file_name, iunit, ierror )

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

  xmin = real ( rxmin - 0.1 * ( rxmax - rxmin ) )
  xmax = real ( rxmax + 0.1 * ( rxmax - rxmin ) )
  ymin = real ( rymin - 0.1 * ( rymax - rymin ) )
  ymax = real ( rymax + 0.1 * ( rymax - rymin ) )

  call ps_file_head ( file_name, xmax, xmin, ymax, ymin )

  call ps_page_head
!
!  Draw the outline of the region.
!
  red = 0.0
  green = 0.0
  blue = 1.0

  call ps_line_rgb ( red, green, blue )

  call ps_moveto ( real ( rxmin ), real ( rymin ) )
  call ps_lineto ( real ( rxmax ), real ( rymin ) )
  call ps_lineto ( real ( rxmax ), real ( rymax ) )
  call ps_lineto ( real ( rxmin ), real ( rymax ) )
  call ps_lineto ( real ( rxmin ), real ( rymin ) )
!
!  Draw a grid in the region.
!
  nx = 11
  ny = 11

  red = 0.5E+00
  green = 0.5E+00
  blue = 0.5E+00

  call ps_grid_cartesian ( real ( rxmin ), real ( rxmax ), nx, &
    real ( rymin ), real ( rymax ), ny )
!
!  Draw the basis function support disks.
!
  red = 0.2
  green = 0.2
  blue = 0.2

  call ps_line_rgb ( red, green, blue )

  red = 0.9
  green = 0.6
  blue = 0.6

  call ps_fill_rgb ( red, green, blue )

  do i = 1, n

    x = center(1,i)
    y = center(2,i)
    r = radius(i)

    if ( filled ) then
      call ps_disk ( x, y, r )
    end if

    call ps_circle ( x, y, r )

  end do

  call ps_page_tail

  call ps_file_tail

  call ps_file_close ( iunit )

  write ( *, * ) ' '
  write ( *, * ) 'BASIS_PLOT'
  write ( *, * ) '  Created a basis function plot file ' // trim ( file_name )

  return
end
subroutine basis_read ( ndim, n, center, radius, input_file )
!
!*******************************************************************************
!
!! BASIS_READ reads basis data from a file.
!
!
!  Modified:
!
!    24 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Output, integer N, the number of points.
!
!    Output, double precision CENTER(NDIM,N), the basis function centers.
!
!    Output, double precision RADIUS(N), the basis function radii.
!
!    Input, character ( len = * ) INPUT_FILE, the name of the file.
!
  implicit none
!
  integer ndim
!
  double precision, dimension ( ndim ) :: c_temp
  double precision, dimension ( ndim, * ) :: center
  integer i
  character ( len = * ) input_file
  integer ios
  integer n
  double precision r_temp
  double precision, dimension ( * ) :: radius
!
  write ( *, * ) ' '
  write ( *, * ) 'BASIS_READ:'
  write ( *, * ) '  Reading basis data file: ', trim ( input_file )

  open ( unit = 12, file = input_file, form = 'formatted', status = 'old' )

  n = 0

  do

    read ( 12, *, iostat = ios ) c_temp(1:ndim), r_temp

    if ( ios /= 0 ) then
      exit
    end if

    n = n + 1
    center(1:ndim,n) = c_temp(1:ndim)
    radius(n) = r_temp

  end do

  close ( unit = 12 )

  write ( *, * ) ' '
  write ( *, * ) 'BASIS_READ:'
  write ( *, * ) '  Number of data points read was ', n

  return
end
subroutine basis_write ( ndim, n, center, radius, output_file )
!
!*******************************************************************************
!
!! BASIS_WRITE writes the basis data to a file.
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, double precision CENTER(NDIM,N), the basis function centers.
!
!    Input, double precision RADIUS(N), the basis function radii.
!
!    Input, character ( len = * ) OUTPUT_FILE, the name of the output file.
!
  implicit none
!
  integer n
  integer ndim
!
  integer i
  character ( len = * ) output_file
  double precision, dimension ( ndim, n ) :: center
  double precision, dimension ( n ) :: radius
!
  write ( *, * ) ' '
  write ( *, * ) 'BASIS_WRITE:'
  write ( *, * ) '  Write basis data file: ', trim ( output_file )
  write ( *, * ) '  Number of data points is ', n

  open ( unit = 12, file = output_file, form = 'formatted', status = 'replace' )

  do i = 1, n
    write ( 12, * ) center(1:ndim,i), radius(i)
  end do

  close ( unit = 12 )

  return
end
subroutine ch_cap ( c )
!
!*******************************************************************************
!
!! CH_CAP capitalizes a single character.
!
!
!  Modified:
!
!    19 July 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character C, the character to capitalize.
!
  implicit none
!
  character c
  integer itemp
!
  itemp = ichar ( c )

  if ( 97 <= itemp .and. itemp <= 122 ) then
    c = char ( itemp - 32 )
  end if

  return
end
subroutine center_plot ( ndim, n, center, radius, a, b )
!
!*******************************************************************************
!
!! CENTER_PLOT plots the center in the region.
!
!
!  Discussion:
!
!    My idea of plotting the basis functions is OK if there are just a
!    few, but quickly becomes useless.  Here, at least, I can plot just
!    the center points, and get an idea of the distribution.
!
!  Modified:
!
!    24 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of basis functions.
!
!    Input, double precision CENTER(NDIM,N), the basis function centers.
!
!    Input, double precision RADIUS(N), the basis function radii.
!
!    Input, double precision A(NDIM), B(NDIM), the minimum and maximum
!    values in each dimension.
!
  implicit none
!
  integer n
  integer ndim
!
  double precision a(ndim)
  double precision b(ndim)
  real blue
  double precision center(ndim,n)
  character ( len = 80 ) file_name
  logical, parameter :: filled = .false.
  real green
  integer i
  integer ierror
  integer iunit
  integer marker_size
  integer nx
  integer ny
  real r
  double precision radius(n)
  real red
  double precision rxmax
  double precision rxmin
  double precision rymax
  double precision rymin
  real x(n)
  real xmax
  real xmin
  real y(n)
  real ymax
  real ymin
!
!  NOT CORRECT FOR NDIM > 2
!
  rxmin = a(1)
  rxmax = b(1)
  rymin = a(2)
  rymax = b(2)

  iunit = 1
  file_name = 'center_plot.ps'

  call ps_file_open ( file_name, iunit, ierror )

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

  xmin = real ( rxmin - 0.1 * ( rxmax - rxmin ) )
  xmax = real ( rxmax + 0.1 * ( rxmax - rxmin ) )
  ymin = real ( rymin - 0.1 * ( rymax - rymin ) )
  ymax = real ( rymax + 0.1 * ( rymax - rymin ) )

  call ps_file_head ( file_name, xmax, xmin, ymax, ymin )

  call ps_page_head
!
!  Draw the outline of the region.
!
  red = 0.0E+00
  green = 0.0E+00
  blue = 1.0E+00

  call ps_line_rgb ( red, green, blue )

  call ps_moveto ( real ( rxmin ), real ( rymin ) )
  call ps_lineto ( real ( rxmax ), real ( rymin ) )
  call ps_lineto ( real ( rxmax ), real ( rymax ) )
  call ps_lineto ( real ( rxmin ), real ( rymax ) )
  call ps_lineto ( real ( rxmin ), real ( rymin ) )
!
!  Draw a grid in the region.
!
  red = 0.2E+00
  green = 0.3E+00
  blue = 0.4E+00

  call ps_line_rgb ( red, green, blue )

  nx = 11
  ny = 11

  call ps_grid_cartesian ( real ( rxmin ), real ( rxmax ), nx, &
    real ( rymin ), real ( rymax ), ny )

  marker_size = 3
  call ps_marker_size ( marker_size )
!
!  Draw the center points.
!
  x(1:n) = real ( center(1,1:n) )
  y(1:n) = real ( center(2,1:n) )

  call ps_mark_circles ( n, x, y )

  call ps_page_tail

  call ps_file_tail

  call ps_file_close ( iunit )

  write ( *, * ) ' '
  write ( *, * ) 'CENTER_PLOT'
  write ( *, * ) '  Created a basis function plot file ' // trim ( file_name )

  return
end
subroutine cvt_make ( alpha, beta, ndim, a, b, maxit, ns, density_function, &
  n, center )
!
!*******************************************************************************
!
!! CVT_MAKE computes centroidal Voronoi tessellation generators.
!
!
!  Discussion:
!
!    The routine is given an initial set of points that define a 
!    Voronoi tessellation.  It computes new points that define a Voronoi
!    tessellation, and which have the property that they are the centroids
!    of their Voronoi regions.
!
!    No stopping criterion is used, for now, except to iterate to MAXIT.
!
!  Modified:
!
!    29 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, double precision ALPHA(2), BETA(2), coefficients used in the
!    generation program.  ALPHA(1) and ALPHA(2) must lie between 0 and 1
!    and sum to 1.  The same is true for BETA.  Common values are
!    ALPHA(1) = BETA(1) = 0,
!    ALPHA(2) = BETA(2) = 1.
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision A(NDIM), B(NDIM), the coordinates of the
!    two extreme corners of the box that defines the region.
!
!    Input, integer MAXIT, the maximum number of iterations.
!
!    Input, integer NS, the average number of sampling points 
!    per generator, on each step of the iteration.
!
!    Input, integer DENSITY_FUNCTION, specifies the density function.
!    1: d(x) = 1.0;
!    2: d(x) = exp ( - 4.0 * ( sum(x(1:n)**2) ) )
!    3: d(x) = exp ( - 3.0 * ( 1.0 - sum(x(1:n)**2) ) )
!
!    Input, integer N, the number of generators.
!
!    Input/output, double precision CENTER(NDIM,N).
!    On input, initial points that generate the Voronoi regions.
!    On output, points that generate the Voronoi regions, which are
!    also the centroids of those regions.
!
!  Local variables:
!
!    Integer UPDATES(N), counts the number of times a generator has been updated.
!
  implicit none
!
  integer n
  integer ndim
!
  double precision a(1:ndim)
  double precision alpha(2)
  double precision b(1:ndim)
  double precision beta(2)
  double precision center(ndim,n)
  double precision center2(ndim,n)
  integer count(n)
  integer density_function
  integer i
  integer ic
  integer iloop
  integer j
  integer k
  integer maxit
  integer ns
  double precision p1(ndim)
  double precision p2(ndim)
  double precision s
  integer updates(n)
  double precision y(ndim)
!
!  UPDATES starts at 1.
!
  updates(1:n) = 1
!
!  Iterate MAXIT times.
!
  do iloop = 1, maxit

    center2(1:ndim,1:n) = 0.0E+00
    count(1:n) = 0

    do j = 1, n * ns
!
!  Generate a random sampling point Y.
!
      call random_generator ( ndim, a, b, density_function, y )
!
!  Find the nearest generator CENTER.  Its index is IC.
!
      call find_closest ( ndim, y, n, center, ic, s )
!
!  Add Y to the averaging data for CENTER(*,IC).
!
      center2(1:ndim,ic) = center2(1:ndim,ic) + y(1:ndim)
      count(ic) = count(ic) + 1

    end do
!
!  For each cell J, average the estimated centroid CENTER with the 
!  averaged hit points CENTER2.  (There are COUNT(J) of these new points).
!
    do j = 1, n

      if ( count(j) /= 0.0 ) then
!
!  This loop and the next could be combined, eliminating the
!  temporary variables P1 and P2.
!
        do k = 1, ndim
          
          p1(k) = ( alpha(1) * updates(j) + beta(1) ) * center(k,j)
          p2(k) = ( alpha(2) * updates(j) + beta(2) ) * center2(k,j) &
            / dble ( count(j) )
        end do

        do k = 1, ndim
          center(k,j) = ( p1(k) + p2(k) ) / dble ( updates(j) + 1 )
        end do

        updates(j) = updates(j) + 1

      end if

    end do

  end do

  return
end
subroutine cvt_read ( ndim, n, center, radius, input_file )
!
!*******************************************************************************
!
!! CVT_READ reads the CVT data from a file.
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Output, double precision CENTER(NDIM,N), the centroidal points.
!
!    Input, character ( len = * ) INPUT_FILE, the name of the file.
!
  implicit none
!
  integer n
  integer ndim
!
  integer i
  character ( len = * ) input_file
  double precision, dimension ( ndim, n ) :: center
  double precision, dimension ( n ) :: radius
!
  write ( *, * ) ' '
  write ( *, *)  'CVT_READ:'
  write ( *, * ) '  Reading Centroidal Voronoi data file: ', trim ( input_file )
  write ( *, * ) '  Number of data points is ', n

  open ( unit = 12, file = input_file, form = 'formatted', status = 'old' )

  do i = 1, n
    read ( 12, * ) center(1:ndim,i), radius(i)
  end do

  close ( unit = 12 )

  return
end
subroutine cvt_write ( ndim, n, center, radius, output_file )
!
!*******************************************************************************
!
!! CVT_WRITE writes the CVT data to a file.
!
!
!  Modified:
!
!    05 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, double precision CENTER(NDIM,N), the centroidal points.
!
!    Input, character ( len = * ) OUTPUT_FILE, the name of the output file.
!
  implicit none
!
  integer n
  integer ndim
!
  integer i
  character ( len = * ) output_file
  double precision, dimension ( ndim, n ) :: center
  double precision, dimension ( n ) :: radius
!
  write ( *, * ) ' '
  write ( *, *)  'CVT_WRITE:'
  write ( *, * ) '  Write Centroidal Voronoi data file: ', trim ( output_file )
  write ( *, * ) '  Number of data points is ', n

  open ( unit = 12, file = output_file, form = 'formatted', status = 'replace' )

  do i = 1, n
    write ( 12, * ) center(1:ndim,i), radius(i)
  end do

  close ( unit = 12 )

  return
end
function density ( ndim, x, density_function )
!
!*******************************************************************************
!
!! DENSITY evaluates the density function.
!
!
!  Discussion:
!
!    The density function is used to control the generation of random points
!    that sample the region.
!
!  Modified:
!
!    19 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, real X(NDIM), the location of the point.
!
!    Input, integer DENSITY_FUNCTION, chooses the density function.
!    1: d(x) = 1.0;
!    2: d(x) = exp ( - 4.0 * ( sum(x(1:n)**2) ) )
!    3: d(x) = exp ( - 3.0 * ( 1.0 - sum(x(1:n)**2) ) )
!
!    Output, real DENSITY, the value of the density function, which
!    should be between 0 and 1 in the region.
!
  implicit none
!
  integer ndim
!
  double precision density
  integer density_function
  double precision, parameter :: pi = &
    3.14159265358979323846264338327950288419716939937510D+00
  double precision x(ndim)
!
!  Here is an odd feature that sets density to zero if the point
!  is more than 1 away from the origin.  This is not what I documented,
!  and should probably be controllable by the user!
!
  if ( sum ( x(1:ndim)**2 ) > 1.0D+00 ) then
    density = 0.0D+00
    return
  end if

  if ( density_function == 1 ) then
    density = 1.0E+00
  else if ( density_function == 2 ) then
    density = exp ( - 4.0E+00 * sum ( x(1:ndim)**2 ) )
  else if ( density_function == 3 ) then
    density = exp ( - 3.0E+00 * ( 1.0E+00 - sum ( x(1:ndim)**2 ) ) )
  else
    density = 1.0D+00
  end if

  return
end
subroutine file_column_count ( file_name, ncolumn )
!
!*******************************************************************************
!
!! FILE_COLUMN_COUNT counts the number of columns in the first line of a file.
!
!
!  Discussion:
!
!    The file is assumed to be a simple text file.
!
!    Each line of the file is presumed to consist of NCOLUMN words, separated
!    by spaces.
!
!    The routine reads the first line and counts the number of words.
!
!  Modified:
!
!    25 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file.
!
!    Output, integer NCOLUMN, the number of columns assumed to be in the file.
!
  implicit none
!
  character ( len = * ) file_name
  integer ios
  integer iunit
  character ( len = 256 ) line
  integer ncolumn
!
!  Open the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'old', form = 'formatted', &
    access = 'sequential', iostat = ios )

  if ( ios /= 0 ) then
    ncolumn = - 1
    write ( *, * ) ' '
    write ( *, * ) 'FILE_LINE_COUNT - Fatal error!'
    write ( *, * ) '  Could not open the file:'
    write ( *, * ) '    ' // trim ( file_name )
    return
  end if
!
!  Read one line.
!
  read ( iunit, '(a)', iostat = ios ) line

  if ( ios /= 0 ) then
    ncolumn = -2
    write ( *, * ) ' '
    write ( *, * ) 'FILE_LINE_COUNT - Fatal error!'
    write ( *, * ) '  The initial line of the file could not be read.'
    write ( *, * ) '    ' // trim ( file_name )
    return
  end if

  close ( unit = iunit )

  call word_count ( line, ncolumn )

  return
end
subroutine file_line_count ( file_name, nline )
!
!*******************************************************************************
!
!! FILE_LINE_COUNT counts the number of lines in a file.
!
!
!  Discussion:
!
!    The file is assumed to be a simple text file.
!
!  Modified:
!
!    21 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file.
!
!    Output, integer NLINE, the number of lines found in the file.
!
  implicit none
!
  character ( len = * ) file_name
  integer ios
  integer iunit
  integer nline
!
  nline = 0
!
!  Open the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'old', form = 'formatted', &
    access = 'sequential', iostat = ios )

  if ( ios /= 0 ) then
    nline = - 1
    write ( *, * ) ' '
    write ( *, * ) 'FILE_LINE_COUNT - Fatal error!'
    write ( *, * ) '  Could not open the file:'
    write ( *, * ) '    ' // trim ( file_name )
    return
  end if
!
!  Count the lines.
!
  do

    read ( iunit, '(a)', iostat = ios )

    if ( ios /= 0 ) then
      exit
    end if

    nline = nline + 1

  end do

  close ( unit = iunit )

  return
end
subroutine find_closest ( ndim, y, n, center, ic, s )
!
!*******************************************************************************
!
!! FIND_CLOSEST finds the center point CENTER(1:NDIM,IC) closest to Y(1:NDIM).
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision Y(NDIM), the point to be checked.
!
!    Input, integer N, the number of center points.
!
!    Input, double precision CENTER(NDIM,N), the center points.
!
!    Output, integer IC, the index of the nearest center point.
!
!    Output, double precision S, the distance.
!
  implicit none
!
  integer n
  integer ndim
!
  double precision center(ndim,n)
  double precision dist_sq
  integer i
  integer ic
  double precision s
  double precision y(ndim)
!
  s = huge ( s )

  do i = 1, n

    dist_sq = sum ( ( center(1:ndim,i) - y(1:ndim) )**2 )

    if ( dist_sq < s ) then
      s = dist_sq
      ic = i
    end if

  end do

  s = sqrt ( s )

  return
end
subroutine find_re ( ndim, pt, r, n, center, np, dist, ord, ierr )
!
!*******************************************************************************
!
!! FIND_RE seems to find the number of center points near a given point.
!
!
!  Modified:
!
!    25 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, real PT(NDIM), the base point.
!
!    Input, real R, the maximum distance from PT to be considered.
!
!    Input, integer N, the number of center points.
!
!    Input, real CENTER(NDIM,N), the center points.
!
!    Output, integer NP, the number of neighboring center points found.
!
!    Output, real DIST(1:NP), lists the Euclidean distance between PT and
!    each of the neighboring center points found.
!
!    Output, integer ORD(1:NP), lists the indices in CENTER of the points found.
!
!    Output, integer IERR, error flag.
!    -1, if no points at all were found.
!    0, if at least one point was found.
!
  implicit none
!
  integer n
  integer ndim
!
  double precision center(ndim,n)
  double precision dist(n)
  double precision dx
  double precision dy
  integer i
  integer ierr
  integer j
  integer np
  integer ord(n)
  double precision pt(ndim)
  double precision r
!
!  Search for a center point CENTER(I) which is no more than R away in X or Y,
!  and which is not identical to PT.  Count the number of such points found,
!  their distance, and indexes.
!
  np = 0

  do i = 1, n
!
!  We do not allow the case where a point in CENTER is exactly equal to PT.
!
    if ( all ( pt(1:ndim) == center(1:ndim,i) ) ) then
      cycle
    end if
!
!  A center point is a suitable neighbor if its L-Infinity distance
!  from PT is no more than R.  But for some reason, we then record
!  its L2 distance...Is this an inconsistency?
!
    if ( all ( abs ( pt(1:ndim) - center(1:ndim,i) ) <= r ) ) then
      np = np + 1
      dist(np) = sqrt ( sum ( ( pt(1:ndim) - center(1:ndim,i) )**2 ) )
      ord(np) = i
    end if

  end do

  if ( np == 0 ) then
    ierr = -1
  else
    ierr = 0
  end if

  return
end
subroutine get_unit ( iunit )
!
!*******************************************************************************
!
!! 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.
!
!  Modified:
!
!    02 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5 and 6).
!
!    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.
!
  implicit none
!
  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 ) 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 halton_make ( ndim, n, points )
!
!*******************************************************************************
!
!! HALTON_MAKE sets up the Halton data points.
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of data points to generate.
!
!    Output, double precision POINTS(NDIM,N), the data points.
!
  implicit none
!
  integer, parameter :: max_prime = 10
  integer n
  integer ndim
!
  integer i
  integer j
  double precision phi
  integer, save, dimension ( max_prime ) :: prime = &
    (/ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 /)
  double precision, parameter :: pi = &
    3.14159265358979323846264338327950288419716939937510D+00
  double precision, dimension ( ndim, n ) :: points
  integer q
  double precision t
!
  if ( ndim > max_prime ) then
    write ( *, * ) ' '
    write ( *, * ) 'HALTON_MAKE - Fatal error!'
    write ( *, * ) '  We did not store enough primes.'
    write ( *, * ) '  NDIM = ', ndim
    write ( *, * ) '  MAX_PRIME = ', max_prime
    stop
  end if

  do i = 1, n
    do j = 1, ndim
      call i_to_halton ( i, prime(j), points(j,i) )
    end do
  end do
!
!  THIS DOESN'T WORK FOR NDIM > 2.
!
!  And besides, I actually want to map this into [A(1:NDIM), B(1:NDIM)]...
!
!  Map a point in the unit square into the unit circle.
!
  do i = 1, n
    t = points(1,i)
    phi = 2.0D+00 * pi * points(2,i)
    points(1,i) = sqrt ( 1.0D+00 - t**2 ) * cos ( phi )
    points(2,i) = sqrt ( 1.0D+00 - t**2 ) * sin ( phi )
  end do

  write ( *, * ) '  Number of Halton data points created was ', n

  return
end
subroutine halton_read ( ndim, n, points, input_file )
!
!*******************************************************************************
!
!! HALTON_READ reads Halton data from a file.
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Output, double precision POINTS(NDIM,N), the Halton points.
!
!    Input, character ( len = * ) INPUT_FILE, the name of the file.
!
  implicit none
!
  integer n
  integer ndim
!
  integer i
  character ( len = * ) input_file
  double precision, dimension ( ndim, n ) :: points
!
  write ( *, * ) ' '
  write ( *, *)  'HALTON_READ:'
  write ( *, * ) '  Reading Halton data file: ', trim ( input_file )
  write ( *, * ) '  Number of data points is ', n

  open ( unit = 12, file = input_file, form = 'formatted', status = 'old' )

  do i = 1, n
    read ( 12, * ) points(1:ndim,i)
  end do

  close ( unit = 12 )

  return
end
subroutine halton_write ( ndim, n, points, output_file )
!
!*******************************************************************************
!
!! HALTON_WRITE writes the Halton data to a file.
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, double precision POINTS(NDIM,N), the Halton points.
!
!    Input, character ( len = * ) OUTPUT_FILE, the name of the output file.
!
  implicit none
!
  integer n
  integer ndim
!
  integer i
  character ( len = * ) output_file
  double precision, dimension ( ndim, n ) :: points
!
  write ( *, * ) ' '
  write ( *, *)  'HALTON_WRITE:'
  write ( *, * ) '  Write Halton data file: ', trim ( output_file )
  write ( *, * ) '  Number of data points is ', n

  open ( unit = 12, file = output_file, form = 'formatted', status = 'replace' )

  do i = 1, n
    write ( 12, * ) points(1:ndim,i)
  end do

  close ( unit = 12 )

  return
end
subroutine i_to_halton ( n, p, r )
!
!*******************************************************************************
!
!! I_TO_HALTON computes an element of a Halton sequence.
!
!
!  Discussion:
!
!    The Halton sequence is often used to generate a "subrandom"
!    sequence of points which have a better covering property
!    than pseudorandom points.
!
!    The Halton sequence generates a sequence of points in [0,1]
!    which (theoretically) never repeats.
!
!    To generate a sequence of points in a 2 dimensional space,
!    it is typical practice to generate a pair of Halton sequences,
!    the first with prime base 2, the second with prime base 3.
!    Similarly, by using the first K primes, a suitable sequence
!    in K-dimensional space can be generated.
!
!    The generation is quite simple.  Given an integer N, the expansion
!    of N in base P is generated.  Then, essentially, the result R
!    is generated by writing a decimal point followed by the digits of
!    the expansion of N, in reverse order.  This decimal value is actually
!    still in base P, so it must be properly interpreted to generate
!    a usable value.
!
!  Example:
!
!    P = 2
!
!    N        N         Halton   Halton
!    decimal  binary    binary   decimal
!    -------  ------    ------   -------
!        1  =     1  =>  .1     = 0.5
!        2  =    10  =>  .01    = 0.25
!        3  =    11  =>  .11    = 0.75
!        4  =   100  =>  .001   = 0.125
!        5  =   101  =>  .101   = 0.625
!        6  =   110  =>  .011   = 0.375
!        7  =   111  =>  .111   = 0.875
!        8  =  1000  =>  .0001  = 0.0625
!
!  Reference:
!
!    J H Halton,
!    Numerische Mathematik,
!    Volume 2, pages 84-90.
!
!  Modified:
!
!    17 December 2000
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer N, the index of the desired element.
!
!    Input, integer P, the Halton base, which should be a prime number.
!
!    Output, double precision R, the N-th element of the Halton sequence
!    for base P.
!
  implicit none
!
  integer digit
  integer n
  integer n2
  integer p
  double precision pp
  double precision r
!
  n2 = abs ( n )
  r = 0.0D+00
  pp = 1.0D+00 / dble ( p )

  do while ( n2 /= 0 )
    digit = mod ( n2, p )
    r = r + digit * pp
    pp = pp / real ( p )
    n2 = n2 / p
  end do

  if ( n < 0 ) then
    r = - r
  end if

  return
end
subroutine radius_make_1 ( ndim, a, b, n, m, density_function, center, radius )
!
!*******************************************************************************
!
!! RADIUS_MAKE_1 uses algorithm 1 to determine a suitable radius for each center.
!
!
!  Modified:
!
!    25 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision A(NDIM), B(NDIM), the coordinates of the
!    two extreme corners of the box that defines the region.
!
!    Input, integer N, the number of Voronoi regions.
!
!    Input, integer M, ?
!
!    Input, integer DENSITY_FUNCTION, specifies the density function.
!    1: d(x) = 1.0;
!    2: d(x) = exp ( - 4.0 * ( sum(x(1:n)**2) ) )
!    3: d(x) = exp ( - 3.0 * ( 1.0 - sum(x(1:n)**2) ) )
!
!    Input, double precision CENTER(NDIM,N), the centers of Voronoi regions.
!
!    Output, double precision RADIUS(N), the radii.
!
  implicit none
!
  integer m
  integer n
  integer ndim
!
  double precision a(ndim)
  double precision asspts(ndim,m*m)
  double precision b(ndim)
  double precision cc
  double precision center(ndim,n)
  double precision dd
  integer density_function
  integer i
  integer ic
  integer mm
  double precision pt(ndim)
  double precision radius(n)
  double precision s
!
!  Is this computation of M*M what we want for NDIM > 2?
!
  mm = m * m

  call random_number ( asspts(1:ndim,1:mm/2) )
!
!  This CAN'T be right for NDIM > 2:
!  And how are we being assured that ASSPTS(2,I) is between A(2) and B(2)???
!
  do i = 1, mm/2
    asspts(1,i) = a(1) + ( b(1) - a(1) ) * asspts(1,i)
    cc = - sqrt ( 1.0D+00 - asspts(1,i)**2 )
    dd = - cc
    asspts(2,i) = cc + ( dd - cc ) * asspts(2,i)
  end do

  do i = mm/2+1, mm
    call random_generator ( ndim, a, b, density_function, asspts(1,i) )
  end do

  radius(1:n) = 0.0

  do i = 1, mm
    pt(1:ndim) = asspts(1:ndim,i)
    call find_closest ( ndim, pt, n, center, ic, s )

    if ( radius(ic) < s ) then
      radius(ic) = s
    end if

  end do

  return
end
subroutine radius_make_2 ( ndim, a, b, n, m, center, radius )
!
!*******************************************************************************
!
!! RADIUS_MAKE_2 uses algorithm 2 to determine a suitable radius for each center.
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision A(NDIM), B(NDIM), the coordinates of the
!    two extreme corners of the box that defines the region.
!
!    Input, integer N, the number of Voronoi regions.
!
!    Input, integer M, ?
!
!    Input, double precision CENTER(NDIM,N), the centers of Voronoi regions.
!
!    Output, double precision RADIUS(N), the radii.
!
  implicit none
!
  integer m
  integer n
  integer ndim
!
  double precision a(ndim)
  double precision b(ndim)
  double precision center(ndim,n)
  double precision dd
  double precision dist(n)
  double precision hx
  double precision hy
  integer i
  integer ierr
  integer j
  integer kk
  integer np
  integer mm
  integer mmn
  integer ord(n)
  double precision pps(ndim,n+m*m)
  double precision pt(ndim)
  double precision r
  double precision radius(n)
!
  if ( m <= 1 ) then
    write ( *, * ) ' '
    write ( *, * ) 'RADIUS_MAKE_2 - Fatal error!'
    write ( *, * ) '  Input parameter M <= 1.'
    write ( *, * ) '  M = ', m
    stop
  end if

  mm = m*m
  mmn = mm + n

  pps(1:ndim,1:n) = center(1:ndim,1:n)
!
!  Select M * M evenly spaced points in [A(1:NDIM),B(1:NDIM)]
!
!  Need to modify this for NDIM > 2
!  This is hard-wired for 2D right now!
!
  do i = 1, m
    do j = 1, m
      pps(1,n+(i-1)*m+j) = a(1) + dble (i-1) * ( b(1) - a(1) ) / dble ( m - 1 )
      pps(2,n+(i-1)*m+j) = a(2) + dble (j-1) * ( b(2) - a(2) ) / dble ( m - 1 )
    end do
  end do
!
!  This is hard-wired for 2D right now!
!
  r = 0.2 * sqrt ( ( b(1) - a(1) ) * ( b(2) - a(2) ) / dble ( n ) )

  radius(1:n) = 0.0

  do i = 1, mm+n

    pt(1:ndim) = pps(1:ndim,i)
    ierr = -1
    dd = 0.0D+00
    kk = 0

    do while ( ierr /= 0 )

      call find_re ( ndim, pt, r, n, center, np, dist, ord, ierr )

      if ( ierr /= 0 ) then
        r = 1.5 * r
      else
        dd = 2.0
        do j = 1, np
          if ( dist(j) < dd ) then
            dd = dist(j)
            kk = ord(j)
          end if
        end do
      end if

    end do

    if ( kk /= 0 ) then
      radius(kk) = max ( radius(kk), dd )
    end if

  end do

  return
end
subroutine random_generator ( ndim, a, b, density_function, z )
!
!*******************************************************************************
!
!! RANDOM_GENERATOR returns a random point Z(1:NDIM) in [A(1:NDIM),B(1:NDIM)].
!
!
!  Discussion:
!
!    The region is an NDIM dimensional box.
!
!    A density function DENSITY(X(1:NDIM) controls whether the points are 
!    uniformly distributed in the region, or are clustered or scattered.
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision A(NDIM), B(NDIM), the lower and upper
!    ranges for the variable.
!
!    Input, integer DENSITY_FUNCTION, specifies the density function.
!    1: density(x) = 1.0
!    2: density(x) = exp ( - 4.0 * ( sum(x(1:n)**2) ) )
!    3: density(x) = exp ( - 3.0 * ( 1.0 - sum(x(1:n)**2) ) )
!
!    Output, double precision Z(NDIM), the random point.
!
  implicit none
!
  integer ndim
!
  double precision a(ndim)
  double precision b(ndim)
  double precision density
  integer density_function
  integer i
  double precision r
  double precision x(ndim)
  double precision z(ndim)
!
  do
!
!  Generate a point at random.
!
    do i = 1, ndim
      call random_number ( r )
      x(i) = ( ( 1.0D+00 - r ) * a(i) + r * b(i) )
    end do
!
!  The density function determines if we should accept or reject the point.
!
    call random_number ( r )

    if ( r < density ( ndim, x, density_function ) ) then
      z(1:ndim) = x(1:ndim)
      exit
    end if

  end do

  return
end
function s_eqi ( s1, s2 )
!
!*******************************************************************************
!
!! S_EQI is a case insensitive comparison of two strings for equality.
!
!
!  Examples:
!
!    S_EQI ( 'Anjana', 'ANJANA' ) is .TRUE.
!
!  Modified:
!
!    14 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S1, S2, the strings to compare.
!
!    Output, logical S_EQI, the result of the comparison.
!
  implicit none
!
  character c1
  character c2
  integer i
  integer len1
  integer len2
  integer lenc
  logical s_eqi
  character ( len = * ) s1
  character ( len = * ) s2
!
  len1 = len ( s1 )
  len2 = len ( s2 )
  lenc = min ( len1, len2 )

  s_eqi = .false.

  do i = 1, lenc

    c1 = s1(i:i)
    c2 = s2(i:i)
    call ch_cap ( c1 )
    call ch_cap ( c2 )

    if ( c1 /= c2 ) then
      return
    end if

  end do

  do i = lenc + 1, len1
    if ( s1(i:i) /= ' ' ) then
      return
    end if
  end do

  do i = lenc + 1, len2
    if ( s2(i:i) /= ' ' ) then
      return
    end if
  end do

  s_eqi = .true.

  return
end
subroutine s_to_i ( s, ival, ierror, last )
!
!*******************************************************************************
!
!! S_TO_I reads an integer value from a string.
!
!
!  Modified:
!
!    28 June 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, a string to be examined.
!
!    Output, integer IVAL, the integer value read from the string.
!    If the string is blank, then IVAL will be returned 0.
!
!    Output, integer IERROR, an error flag.
!    0, no error.
!    1, an error occurred.
!
!    Output, integer LAST, the last character of S used to make IVAL.
!
  implicit none
!
  character c
  integer i
  integer ierror
  integer isgn
  integer istate
  integer ival
  integer last
  character ( len = * ) s
!
  ierror = 0
  istate = 0
  isgn = 1
  ival = 0

  do i = 1, len_trim ( s )

    c = s(i:i)
!
!  Haven't read anything.
!
    if ( istate == 0 ) then

      if ( c == ' ' ) then

      else if ( c == '-' ) then
        istate = 1
        isgn = -1
      else if ( c == '+' ) then
        istate = 1
        isgn = + 1
      else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
        istate = 2
        ival = ichar ( c ) - ichar ( '0' )
      else
        ierror = 1
        return
      end if
!
!  Have read the sign, expecting digits.
!
    else if ( istate == 1 ) then

      if ( c == ' ' ) then

      else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
        istate = 2
        ival = ichar ( c ) - ichar ( '0' )
      else
        ierror = 1
        return
      end if
!
!  Have read at least one digit, expecting more.
!
    else if ( istate == 2 ) then

      if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
        ival = 10 * ival + ichar ( c ) - ichar ( '0' )
      else
        ival = isgn * ival
        last = i - 1
        return
      end if

    end if

  end do
!
!  If we read all the characters in the string, see if we're OK.
!
  if ( istate == 2 ) then
    ival = isgn * ival
    last = len_trim ( s )
  else
    ierror = 1
    last = 0
  end if

  return
end
subroutine set_random_seed ( myrank )
!
!*******************************************************************************
!
!! SET_RANDOM_SEED initializes the FORTRAN 90 random number generator.
!
!
!  Modified:
!
!    25 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer MYRANK, an identifier for each processor.  (This
!    is only used in multi-processor applications).
!
  implicit none
!
  character ( len = 10 ) big_ben(3)
  integer date_time(8)
  integer i
  integer k
  integer myrank
  integer, allocatable :: seed(:)
!
!  Initialize the random seed routine.
!
  call random_seed
!
!  Request the size of a typical seed.
!  (It's probably just 1!)
!
  call random_seed ( size = k )
!
!  Set up space for a seed.
!
  allocate ( seed(k) )
!
!  Get the date and time.
!
  call date_and_time ( big_ben(1), big_ben(2), big_ben(3), date_time )
!
!  Make up a "random" value based on date and time information.
!
  do i = 1, k
    seed(i) = date_time(8-mod(i-1,8)) + i * ( myrank + 1 ) * 100
  end do
!
!  Send this random value back to the RANDOM_SEED routine, to be
!  used as the seed of the random number generator.
!
  call random_seed ( put = seed(1:k) )

  return
end
subroutine uniform_make ( ndim, a, b, n, density_function, points )
!
!*******************************************************************************
!
!! UNIFORM_MAKE sets up the uniform random data points.
!
!
!  Modified:
!
!    25 January 2001
!
!  Author:
!
!    Lili Ju
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision A(NDIM), B(NDIM), the coordinates of the
!    two extreme corners of the box that defines the region.
!
!    Input, integer N, the number of data points to generate.
!
!    Input, integer DENSITY_FUNCTION, specifies the density function.
!    1: d(x) = 1.0
!    2: d(x) = exp ( - 4.0 * ( sum(x(1:n)**2) ) )
!    3: d(x) = exp ( - 3.0 * ( 1.0 - sum(x(1:n)**2) ) )
!
!    Output, double precision POINTS(NDIM,N), the data points.
!
  implicit none
!
  integer n
  integer ndim
!
  double precision a(ndim)
  double precision b(ndim)
  integer density_function
  integer i
  double precision points(ndim,n)
!
  do i = 1, n
    call random_generator ( ndim, a, b, density_function, points(1,i) )
  end do

  write ( *, * ) '  Number of Uniform random data points created was ', n

  return
end
subroutine uniform_read ( ndim, n, points, input_file )
!
!*******************************************************************************
!
!! UNIFORM_READ reads uniform random data from a file.
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Output, double precision POINTS(NDIM,N), the uniform random points.
!
!    Input, character ( len = * ) INPUT_FILE, the name of the file.
!
  implicit none
!
  integer n
  integer ndim
!
  integer i
  character ( len = * ) input_file
  double precision, dimension ( ndim, n ) :: points
!
  write ( *, * ) ' '
  write ( *, *)  'UNIFORM_READ:'
  write ( *, * ) '  Reading Uniform Random data file: ', trim ( input_file )
  write ( *, * ) '  Number of data points is ', n

  open ( unit = 12, file = input_file, form = 'formatted', status = 'old' )

  do i = 1, n
    read ( 12, * ) points(1:ndim,i)
  end do

  close ( unit = 12 )

  return
end
subroutine uniform_write ( ndim, n, points, output_file )
!
!*******************************************************************************
!
!! UNIFORM_WRITE writes the uniform random data to a file.
!
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, double precision POINTS(NDIM,N), the uniform random points.
!
!    Input, character ( len = * ) OUTPUT_FILE, the name of the output file.
!
  implicit none
!
  integer n
  integer ndim
!
  integer i
  character ( len = * ) output_file
  double precision, dimension ( ndim, n ) :: points
!
  write ( *, * ) ' '
  write ( *, *)  'UNIFORM_WRITE:'
  write ( *, * ) '  Write Uniform Random data file: ', trim ( output_file )
  write ( *, * ) '  Number of data points is ', n

  open ( unit = 12, file = output_file, form = 'formatted', status = 'replace' )

  do i = 1, n
    write ( 12, * ) points(1:ndim,i)
  end do

  close ( unit = 12 )

  return
end
subroutine word_count ( s, nword )
!
!*******************************************************************************
!
!! WORD_COUNT counts the number of "words" in a string.
!
!
!  Modified:
!
!    14 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, the string to be examined.
!
!    Output, integer NWORD, the number of "words" in the string.
!    Words are presumed to be separated by one or more blanks.
!
  implicit none
!
  logical blank
  integer i
  integer lens
  integer nword
  character ( len = * ) s
!
  nword = 0
  lens = len ( s )

  if ( lens <= 0 ) then
    return
  end if

  blank = .true.

  do i = 1, lens

    if ( s(i:i) == ' ' ) then
      blank = .true.
    else if ( blank ) then
      nword = nword + 1
      blank = .false.
    end if

  end do

  return
end