subroutine create_brownian ( m, n, x )
!
!*******************************************************************************
!
!! CREATE_BROWNIAN creates a random sequence of Brownian points.
!
!
!  Discussion:
!
!    A starting point is generated at the origin.  The next point
!    is generated at a uniformaly random angle and a (0,1) normally
!    distributed distance from the previous point.
!
!    It is up to the user to rescale the data, if desired.
!
!  Modified:
!
!    12 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the dimension of the space.
!
!    Input, integer N, the number of points to create and write.
!
!    Output, real X(M,N), the random points.
!
  implicit none
!
  integer m
  integer n
!
  real direction(m)
  integer j
  real r
  real x(m,n)
!
!  Initial point.
!
  x(1:m,1) = 0.0E+00
!
!  Generate angles and steps.
!
  do j = 2, n

    call normal_01_sample ( r )
    r = abs ( r )

    call direction_random_nd ( m, direction )

    x(1:m,j) = x(1:m,j-1) + r * direction(1:m)

  end do

  return
end
subroutine create_normal ( m, n, x )
!
!*******************************************************************************
!
!! CREATE_NORMAL creates N normally distributed points in M space.
!
!
!  Modified:
!
!    12 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the dimension of the space.
!
!    Input, integer N, the number of points to create and write.
!
!    Output, real X(M,N), the random points.
!
  implicit none
!
  integer m
  integer n
!
  real x(m,n)
!
  call normal_01_vector ( m*n, x(1:m,1:n) )

  return
end
subroutine create_uniform_annulus01 ( m, n, x, r )
!
!*******************************************************************************
!
!! CREATE_UNIFORM_ANNULUS01 creates a random set of points in a unit annulus.
!
!
!  Discussion:
!
!    The annulus is the area between two concentric spheres.  The
!    inner sphere has radius R, and the outer sphere has radius 1.
!
!  Modified:
!
!    10 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the dimension of the space.
!
!    Input, integer N, the number of points to create and write.
!
!    Output, real X(M,N), the random points.
!
!    Input, real R, the inner radius of the annulus, which must
!    be less than 1.
!
  implicit none
!
  integer m
  integer n
!
  integer j
  real r
  real x(m,n)
!
  if ( r**2 >= 1.0E+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CREATE_UNIFORM_ANNULUS_01 - Fatal error!'
    write ( *, '(a)' ) '  R**2 >= 1.'
    return
  end if
!
!  Generate the points using rejection.
!
  do j = 1, n

    do

      call random_number ( harvest = x(1:m,j) )

      x(1:m,j) = 2.0E+00 * x(1:m,j) - 1.0E+00

      if ( r**2 <= sum ( x(1:m,j)**2 ) .and. &
        sum ( x(1:m,j)**2 ) <= 1.0E+00 ) then
        exit
      end if

    end do

  end do

  return
end
subroutine create_uniform_sphere01 ( m, n, x )
!
!*******************************************************************************
!
!! CREATE_UNIFORM_SPHERE_01 creates a random set of points in the unit sphere.
!
!
!  Discussion:
!
!    The sphere has radius 1.
!
!  Modified:
!
!    10 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the dimension of the space.
!
!    Input, integer N, the number of points to create and write.
!
!    Output, real X(M,N), the random points.
!
  implicit none
!
  integer m
  integer n
!
  integer j
  real x(m,n)
  real xx
  real yy
!
  do j = 1, n

    do

      call random_number ( harvest = x(1:m,j) )

      x(1:m,j) = 2.0E+00 * x(1:m,j) - 1.0E+00

      if ( sum ( x(1:m,j)**2 ) <= 1.0E+00 ) then
        exit
      end if

    end do

  end do

  return
end
subroutine create_uniform_block01 ( m, n, x )
!
!*******************************************************************************
!
!! CREATE_UNIFORM_BLOCK01 creates a random set of points in the unit block.
!
!
!  Discussion:
!
!    The block is defined as points all of whose components are between
!    0 and 1.
!
!  Modified:
!
!    10 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the dimension of the space.
!
!    Input, integer N, the number of points to create and write.
!
!    Output, real X(M,N), the random points.
!
  implicit none
!
  integer m
  integer n
!
  real x(m,n)
!
  call random_number ( harvest = x(1:m,1:n) )

  return
end
subroutine direction_random_nd ( m, w )
!
!*******************************************************************************
!
!! DIRECTION_RANDOM_ND generates a random direction vector in ND.
!
!
!  Modified:
!
!    17 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the dimension of the space.
!
!    Output, real W(M), a random direction vector, with unit norm.
!
  implicit none
!
  integer m
!
  integer i
  real norm
  real w(m)
!
!  Get N values from a standard normal distribution.
!
  call normal_01_vector ( m, w )
!
!  Compute the length of the vector.
!
  norm = sqrt ( sum ( w(1:m)**2 ) )
!
!  Normalize the vector.
!
  w(1:m) = w(1:m) / norm

  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 normal_01_vector ( n, x )
!
!*******************************************************************************
!
!! NORMAL_01_VECTOR samples the standard normal probability distribution.
!
!
!  Discussion:
!
!    The standard normal probability distribution function (PDF) has
!    mean 0 and standard deviation 1.
!
!    This routine can generate a vector of values on one call.  It
!    has the feature that it should provide the same results
!    in the same order no matter how we break up the task.
!
!    Before calling this routine, the user may call RANDOM_SEED
!    in order to set the seed of the random number generator.
!
!  Method:
!
!    The Box-Muller method is used, which is efficient, but
!    generates an even number of values each time.  On any call
!    to this routine, an even number of new values are generated.
!    Depending on the situation, one value may be left over.
!    In that case, it is saved for the next call.
!
!  Modified:
!
!    03 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of values desired.  If N is negative,
!    then the code will flush its internal memory; in particular,
!    if there is a saved value to be used on the next call, it is
!    instead discarded.  This is useful if the user has reset the
!    random number seed, for instance.
!
!    Output, real X(N), a sample of the standard normal PDF.
!
  implicit none
!
  integer n
!
  integer m
  real, parameter :: pi = &
    3.14159265358979323846264338327950288419716939937510E+00
  real r(n+1)
  integer, save :: made = 0
  integer, save :: saved = 0
  real x(n)
  integer xhi
  integer xlo
  real, save :: y = 0.0E+00
!
!  I'd like to allow the user to reset the random number seed.
!  But this won't work properly if we have a saved value Y.
!  I'm making a crock option that allows the user to signal
!  explicitly that any internal memory should be flushed,
!  by passing in a negative value for N.
!
  if ( n < 0 ) then
    n = made
    made = 0
    saved = 0
    y = 0.0E+00
    return
  else if ( n == 0 ) then
    return
  end if
!
!  Record the range of X we need to fill in.
!
  xlo = 1
  xhi = n
!
!  Use up the old value, if we have it.
!
  if ( saved == 1 ) then
    x(1) = y
    saved = 0
    xlo = 2
  end if
!
!  If we don't need any more values, return.
!
  if ( xhi - xlo + 1 == 0 ) then

    return
!
!  If we need just one new value, do that here to avoid null arrays.
!
  else if ( xhi - xlo + 1 == 1 ) then

    call random_number ( harvest = r(1:2) )

    x(xhi) = sqrt ( -2.0E+00 * log ( r(1) ) ) * cos ( 2.0E+00 * pi * r(2) )
    y =      sqrt ( -2.0E+00 * log ( r(1) ) ) * sin ( 2.0E+00 * pi * r(2) )

    saved = 1

    made = made + 2
!
!  If we require an even number of values, that's easy.
!
  else if ( mod ( xhi-xlo+1, 2 ) == 0 ) then

    m = ( xhi-xlo+1 ) / 2

    call random_number ( harvest = r(1:2*m) )

    x(xlo:xhi-1:2) = sqrt ( -2.0E+00 * log ( r(1:2*m-1:2) ) ) &
      * cos ( 2.0E+00 * pi * r(2:2*m:2) )

    x(xlo+1:xhi:2) = sqrt ( -2.0E+00 * log ( r(1:2*m-1:2) ) ) &
      * sin ( 2.0E+00 * pi * r(2:2*m:2) )

    made = made + xhi - xlo + 1
!
!  If we require an odd number of values, we generate an even number,
!  and handle the last pair specially, storing one in X(N), and
!  saving the other for later.
!
  else

    xhi = xhi - 1

    m = ( xhi-xlo+1 ) / 2 + 1

    call random_number ( harvest = r(1:2*m) )

    x(xlo:xhi-1:2) = sqrt ( -2.0E+00 * log ( r(1:2*m-3:2) ) ) &
      * cos ( 2.0E+00 * pi * r(2:2*m-2:2) )

    x(xlo+1:xhi:2) = sqrt ( -2.0E+00 * log ( r(1:2*m-3:2) ) ) &
      * sin ( 2.0E+00 * pi * r(2:2*m-2:2) )

    x(n) = sqrt ( -2.0E+00 * log ( r(2*m-1) ) ) &
      * cos ( 2.0E+00 * pi * r(2*m) )

    y = sqrt ( -2.0E+00 * log ( r(2*m-1) ) ) &
      * sin ( 2.0E+00 * pi * r(2*m) )

    saved = 1

    made = made + xhi - xlo + 2

  end if

  return
end
subroutine normal_01_sample ( x )
!
!*******************************************************************************
!
!! NORMAL_01_SAMPLE samples the standard normal probability distribution.
!
!
!  Discussion:
!
!    The standard normal probability distribution function (PDF) has
!    mean 0 and standard deviation 1.
!
!  Method:
!
!    The Box-Muller method is used, which is efficient, but
!    generates two values at a time.
!
!  Modified:
!
!    02 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real X, a sample of the standard normal PDF.
!
  implicit none
!
  real, parameter :: pi = &
    3.14159265358979323846264338327950288419716939937510E+00
  real r1
  real r2
  integer, save :: used = -1
  real x
  real, save :: y = 0.0E+00
!
  if ( used == -1 ) then
    call random_seed ( )
    used = 0
  end if
!
!  If we've used an even number of values so far, generate two more,
!  return one and save one.
!
  if ( mod ( used, 2 ) == 0 ) then

    do

      call random_number ( harvest = r1 )

      if ( r1 /= 0.0E+00 ) then
        exit
      end if

    end do

    call random_number ( harvest = r2 )

    x = sqrt ( -2.0E+00 * log ( r1 ) ) * cos ( 2.0E+00 * pi * r2 )
    y = sqrt ( -2.0E+00 * log ( r1 ) ) * sin ( 2.0E+00 * pi * r2 )
!
!  Otherwise, return the second, saved, value.
!
  else

    x = y

  end if

  used = used + 1

  return
end
subroutine rescale_sphere01 ( m, n, x )
!
!*******************************************************************************
!
!! RESCALE_CIRCLE01 translates and rescales data to fit in the unit sphere.
!
!
!  Modified:
!
!    11 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the dimension of the space.
!
!    Input, integer N, the number of points to create and write.
!
!    Input/output, real X(M,N), the data to be modified.
!
  implicit none
!
  integer m
  integer n
!
  integer i
  real r(m)
  real scale
  real x(m,n)
  real xave(m)
!
!  Determine the center.
!
  do i = 1, m
    xave(i) = sum ( x(i,1:n) ) / real ( n )
  end do
!
!  Determine the maximum distance from the center.
!
  do i = 1, m
    r(1:m) = sum ( sqrt ( ( x(i,1:n) - xave(i) )**2 ) )
  end do

  scale = maxval ( r(1:m) )

  if ( scale > 0.0E+00 ) then
    do i = 1, m
      x(i,1:n) = 0.5E+00 + ( x(i,1:n) - xave(i) ) / scale
    end do
  else
    x(1:m,1:n) = 0.5E+00
  end if

  return
end
subroutine rescale_block01 ( m, n, x )
!
!*******************************************************************************
!
!! RESCALE_BLOCK01 translates and rescales data to fit in the unit block.
!
!
!  Modified:
!
!    11 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the dimension of the space.
!
!    Input, integer N, the number of points to create and write.
!
!    Input/output, real X(M,N), the data to be modified.
!
  implicit none
!
  integer m
  integer n
!
  integer i
  real x(m,n)
  real xmax(m)
  real xmin(m)
!
  do i = 1, m
    xmin(i) = minval ( x(i,1:n) )
    xmax(i) = maxval ( x(i,1:n) )
  end do

  do i = 1, m
    if ( xmax(i) - xmin(i) > 0.0E+00 ) then
      x(i,1:n) = ( x(i,1:n) - xmin(i) ) / ( xmax(i) - xmin(i) )
    else
      x(i,1:n) = 0.5E+00
    end if
  end do

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

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

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

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

  return
end
subroutine write_close ( iunit )
!
!*******************************************************************************
!
!! WRITE_CLOSE closes a file that was written.
!
!
!  Modified:
!
!    10 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IUNIT, the unit number.
!
  implicit none
!
  integer iunit
!
  close ( unit = iunit )

  return
end
subroutine write_data ( iunit, m, n, x )
!
!*******************************************************************************
!
!! WRITE_DATA writes out a vector of X, Y data.
!
!
!  Modified:
!
!    10 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IUNIT, the unit number.
!
!    Input, integer N, the number of points to write.
!
!    Input, real X(M,N), the data to write.  Each record will
!    contain one set of data, X(1:M,J).
!
  implicit none
!
  integer m
  integer n
!
  integer iunit
  integer j
  real x(m,n)
!
  do j = 1, n
    write ( unit = iunit, '(10f14.8)' ) x(1:m,j)
  end do

  return
end
subroutine write_open ( iunit, file_name )
!
!*******************************************************************************
!
!! WRITE_OPEN opens a file to be written.
!
!
!  Modified:
!
!    10 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IUNIT, the unit number.
!
!    Input, character ( len = * ) FILE_NAME, the name of the file.
!
  implicit none
!
  character ( len = * ) file_name
  integer ios
  integer iunit
!
  call get_unit ( iunit )

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

  if ( ios /= 0 ) then
  end if

  return
end