#! /usr/bin/env python3
#
def polygon_centroid ( n, v ):

#*****************************************************************************80
#
## polygon_centroid() computes the centroid of a polygon.
#
#  Discussion:
#
#    Denoting the centroid coordinates by CENTROID, then
#
#      CENTROID(1) = Integral ( Polygon interior ) x dx dy / Area ( Polygon )
#      CENTROID(2) = Integral ( Polygon interior ) y dx dy / Area ( Polygon ).
#
#    Green's theorem states that for continuously differentiable functions
#    M(x,y) and N(x,y),
#
#      Integral ( Polygon boundary ) ( M dx + N dy ) =
#      Integral ( Polygon interior ) ( dN/dx - dM/dy ) dx dy.
#
#    Using M(x,y) = 0 and N(x,y) = x^2/2, we get:
#
#      CENTROID(1) = 0.5 * Integral ( Polygon boundary ) x^2 dy 
#                  / Area ( Polygon ),
#
#    which becomes
#
#      CENTROID(1) = 1/6 sum ( 1 <= I <= N )
#        ( X(I+1) + X(I) ) * ( X(I) * Y(I+1) - X(I+1) * Y(I))
#        / Area ( Polygon )
#
#    where, when I = N, the index "I+1" is replaced by 1.
#
#    A similar calculation gives us a formula for CENTROID(2).
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    22 May 2005
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Gerard Bashein and Paul Detmer,
#    Centroid of a Polygon,
#    Graphics Gems IV, 
#    edited by Paul Heckbert,
#    AP Professional, 1994.
#
#  Input:
#
#    integer N, the number of sides of the polygon.
#
#    real V(N,2), the coordinates of the vertices.
#
#  Output:
#
#    real CENTROID(2), the coordinates of the centroid.
#
  import numpy as np

  area = 0.0
  centroid = np.zeros ( 2 )

  for i in range ( 0, n ):

    if ( i < n - 1 ):
      ip1 = i + 1
    else:
      ip1 = 0

    temp = ( v[i,0] * v[ip1,1] - v[ip1,0] * v[i,1] )

    area = area + temp

    centroid = centroid + ( v[ip1,:] + v[i,:] ) * temp

  area = area / 2.0

  if ( area == 0.0 ):
    centroid = v[0,:]
  else:
    centroid = centroid / ( 6.0 * area )

  return centroid

def random_data_test ( ):

#*****************************************************************************80
#
## random_data_test() tests random_data().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    25 April 2022
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'random_data_test():' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  random_data() samples points in various geometric objects.' )

  uniform_in_annulus_test ( )
  uniform_in_annulus_sector_test ( )
  uniform_in_circle_test ( )
  uniform_in_ellipse_test ( )
  uniform_in_ellipsoid_test ( )
  uniform_in_hexagon_test ( )
  uniform_in_hypercube_test ( )
  uniform_in_hypersphere_test ( )
  uniform_in_parallelogram_test ( )
  uniform_in_polygon_test ( )
  uniform_in_sector_test ( )
  uniform_in_simplex_test ( )
  uniform_in_tetrahedron_test ( )
  uniform_in_triangle_test ( )

  uniform_on_circle_test ( )
  uniform_on_ellipse_test ( )
  uniform_on_ellipsoid_test ( )
  uniform_on_hemisphere_phong_test ( )
  uniform_on_hypercube_test ( )
  uniform_on_hypersphere_test ( )
  uniform_on_simplex_test ( )
  uniform_on_sphere_patch_test ( )
  uniform_on_sphere_triangle_test ( )
  uniform_on_triangle_test ( )

  uniform_walk_test ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'random_data_test():' )
  print ( '  Normal end of execution.' )

def sphere_triangle_angles_to_area ( r, a, b, c ):

#*****************************************************************************80
#
## sphere_triangle_angles_to_area() computes the area of a spherical triangle.
#
#  Discussion:
#
#    A sphere centered at 0 in 3D satisfies the equation:
#
#      X*X + Y*Y + Z*Z = R*R
#
#    A spherical triangle is specified by three points on the surface
#    of the sphere.
#
#    The area formula is known as Girard's formula.
#
#    The area of a spherical triangle is:
#
#      AREA = ( A + B + C - PI ) * R*R
#
#    where A, B and C are the (surface) angles of the triangle.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 February 2005
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real R, the radius of the sphere.
#
#    real A, B, C, the angles of the triangle.
#
#  Output:
#
#    real AREA, the area of the spherical triangle.
#
  import numpy as np
#
#  Apply Girard's formula.
#
  area = r * r * ( a + b + c - np.pi )

  return area

def sphere_triangle_sides_to_angles ( r, aside, bside, cside ):

#*****************************************************************************80
#
## sphere_triangle_sides_to_angles() computes spherical triangle angles in 3D.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    19 July 2018
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real R, the radius of the sphere.
#
#    real ASIDE, BSIDE, CSIDE, the (geodesic) length of the 
#    sides of the triangle.
#
#  Output:
#
#    real A, B, C, the spherical angles of the triangle.
#    Angle A is opposite the side of length ASIDE, and so on.
#
  import numpy as np

  asu = aside / r
  bsu = bside / r
  csu = cside / r
  ssu = ( asu + bsu + csu ) / 2.0

  tan_a2 = np.sqrt ( ( np.sin ( ssu - bsu ) * np.sin ( ssu - csu ) ) / \
                     ( np.sin ( ssu )       * np.sin ( ssu - asu ) ) )

  a = 2.0 * np.arctan ( tan_a2 )

  tan_b2 = np.sqrt ( ( np.sin ( ssu - asu ) * np.sin ( ssu - csu ) ) / \
                     ( np.sin ( ssu )       * np.sin ( ssu - bsu ) ) )

  b = 2.0 * np.arctan ( tan_b2 )

  tan_c2 = np.sqrt ( ( np.sin ( ssu - asu ) * np.sin ( ssu - bsu ) ) / \
                     ( np.sin ( ssu )       * np.sin ( ssu - csu ) ) )

  c = 2.0 * np.arctan ( tan_c2 )

  return a, b, c

def sphere_triangle_vertices_to_sides ( r, v1, v2, v3 ):

#*****************************************************************************80
#
## sphere_triangle_vertices_to_sides() computes spherical triangle sides.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    25 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real R, the radius of the sphere.
#
#    real V1(3), V2(3), V3(3), the vertices of the spherical triangle.
#
#  Output:
#
#    real ASIDE, BSIDE, CSIDE, the (geodesic) length of the sides
#    of the triangle.
#
  import numpy as np

  aside = r * np.arccos ( np.dot ( v2, v3 ) / r ** 2 )
  bside = r * np.arccos ( np.dot ( v3, v1 ) / r ** 2 )
  cside = r * np.arccos ( np.dot ( v1, v2 ) / r ** 2 )

  return aside, bside, cside

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

def triangle_area ( t ):

#*****************************************************************************80
#
## triangle_area() computes the area of a triangle.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real T(3,2), the triangle vertices.
#
#  Output:
#
#    real AREA, the absolute area of the triangle.
#
  import numpy as np

  area = 0.5 * np.abs ( \
      t[0,0] * ( t[1,1] - t[2,1] ) \
    + t[1,0] * ( t[2,1] - t[0,1] ) \
    + t[2,0] * ( t[0,1] - t[1,1] ) )

  return area

def uniform_in_annulus ( n, pc, r1, r2 ):

#*****************************************************************************80
#
## uniform_in_annulus() samples a circular annulus.
#
#  Discussion:
#
#    A circular annulus with center PC, inner radius R1 and
#    outer radius R2, is the set of points P so that
#
#      R1^2 <= (P(1)-PC(1))^2 + (P(2)-PC(2))^2 <= R2^2
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Peter Shirley,
#    Nonuniform Random Point Sets Via Warping,
#    Graphics Gems, Volume III,
#    edited by David Kirk,
#    AP Professional, 1992, 
#    ISBN: 0122861663,
#    LC: T385.G6973.
#
#  Input:
#
#    integer N, the number of points to generate.
#
#    real PC(2), the center.
#
#    real R1, R2, the inner and outer radii.
#
#  Output:
#
#    real P(N,2), sample points.
#
  import numpy as np

  t = np.random.rand ( n )
  t = 2.0 * np.pi * t

  r = np.random.rand ( n )

  r = np.sqrt ( ( 1.0 - r ) * r1 * r1 \
       +                r   * r2 * r2 )

  p = np.zeros ( [ n, 2 ] )
  p[:,0] = pc[0] + r * np.cos ( t )
  p[:,1] = pc[1] + r * np.sin ( t )

  return p

def uniform_in_annulus_test ( ):

#*****************************************************************************80
#
## uniform_in_annulus_test() tests uniform_in_annulus().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 400
  pc = np.array ( [ 10.0, 5.0 ] )
  r1 = 1.0
  r2 = 3.0

  print ( '' )
  print ( 'uniform_in_annulus_test():' )
  print ( '  uniform_in_annulus() returns uniform random points' )
  print ( '  from inside a circular annulus.' )
  print ( '' )
  print ( '  Number of points N =          ', n )
  print ( '  Inner radius =                ', r1 )
  print ( '  Outer radius =                ', r2 )
  print ( '  Annulus center =              ', pc )

  x = uniform_in_annulus ( n, pc, r1, r2 )

  nc = 101
  t = np.linspace ( 0.0, 2.0 * np.pi, nc )
  x1 = np.zeros ( [ nc, 2 ] )
  x1[:,0] = pc[0] + r1 * np.cos ( t )
  x1[:,1] = pc[1] + r1 * np.sin ( t )
  x2 = np.zeros ( [ nc, 2 ] )
  x2[:,0] = pc[0] + r2 * np.cos ( t )
  x2[:,1] = pc[1] + r2 * np.sin ( t )

  plt.clf ( )
  plt.plot ( x1[:,0], x1[:,1], 'r', linewidth = 3 )
  plt.plot ( x2[:,0], x2[:,1], 'r', linewidth = 3 )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_annulus()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_in_annulus.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_annulus_sector ( n, pc, r1, r2, theta1, theta2 ):

#*****************************************************************************80
#
## uniform_in_annulus_sector() samples an annular sector in 2D.
#
#  Discussion:
#
#    An annular sector with center PC, inner radius R1 and
#    outer radius R2, and angles THETA1, THETA2, is the set of points
#    P so that
#
#      R1^2 <= (P(1)-PC(1))^2 + (P(2)-PC(2))^2 <= R2^2
#
#    and
#
#      THETA1 <= THETA ( P - PC ) <= THETA2
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Peter Shirley,
#    Nonuniform Random Point Sets Via Warping,
#    Graphics Gems, Volume III,
#    edited by David Kirk,
#    AP Professional, 1992, 
#    ISBN: 0122861663,
#    LC: T385.G6973.
#
#  Input:
#
#    integer N, the number of points to generate.
#
#    real PC(2), the center.
#
#    real R1, R2, the inner and outer radii.
#
#    real THETA1, THETA2, the angles.
#
#  Output:
#
#    real P(N,2), sample points.
#
  import numpy as np

  u = np.random.rand ( n )

  theta = ( 1.0 - u ) * theta1 \
     +            u   * theta2

  v = np.random.rand ( n )

  r = np.sqrt ( ( 1.0 - v ) * r1 * r1 \
         +              v   * r2 * r2 )

  p = np.zeros ( [ n, 2 ] )
  p[:,0] = pc[0] + r * np.cos ( theta )
  p[:,1] = pc[1] + r * np.sin ( theta )

  return p

def uniform_in_annulus_sector_test ( ):

#*****************************************************************************80
#
## uniform_in_annulus_sector_test() tests uniform_in_annulus_sector().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 400
  pc = np.array ( [ 10.0, 5.0 ] )
  r1 = 1.0
  r2 = 3.0
  theta1 = 0.0
  theta2 = np.pi / 2.0

  print ( '' )
  print ( 'uniform_in_annulus_sector_test():' )
  print ( '  uniform_in_annulus_sector() returns uniform random points' )
  print ( '  from inside a sector of a circular annulus.' )
  print ( '' )
  print ( '  Number of points N =          ', n )
  print ( '  Inner radius =                ', r1 )
  print ( '  Outer radius =                ', r2 )
  print ( '  Annulus center =              ', pc )
  print ( '  theta1 =                      ', theta1 )
  print ( '  theta2 =                      ', theta2 )

  x = uniform_in_annulus_sector ( n, pc, r1, r2, theta1, theta2 )
#
#  Plot.
#
  nc = 101
  t = np.linspace ( theta1, theta2, nc )
  x1 = np.zeros ( [ nc, 2 ] )
  x1[:,0] = pc[0] + r1 * np.cos ( t )
  x1[:,1] = pc[1] + r1 * np.sin ( t )
  x2 = np.zeros ( [ nc, 2 ] )
  x2[:,0] = pc[0] + r2 * np.cos ( t )
  x2[:,1] = pc[1] + r2 * np.sin ( t )

  plt.clf ( )
  plt.plot ( x1[:,0], x1[:,1], 'r', linewidth = 3 )
  plt.plot ( x2[:,0], x2[:,1], 'r', linewidth = 3 )
  plt.plot ( [x1[0,0], x2[0,0]], [x1[0,1], x2[0,1]], 'r', linewidth = 3 )
  plt.plot ( [x1[-1,0], x2[-1,0]], [x1[-1,1], x2[-1,1]], 'r', linewidth = 3 )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_annulus_sector()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_in_annulus_sector.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_circle ( n ):

#*****************************************************************************80
#
## uniform_in_circle() maps uniform points into the unit circle.
#
#  Discussion:
#
#    The unit circle has center at the origin, and radius 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    03 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of points.
#
#  Output:
#
#    real X(N,2), the points.
#
  import numpy as np

  r = np.random.rand ( n )
  r = np.sqrt ( r )

  t = np.random.rand ( n )
  t = 2.0 * np.pi * t

  x = np.zeros ( [ n, 2 ] )
  x[:,0] = r * np.cos ( t )
  x[:,1] = r * np.sin ( t )

  return x

def uniform_in_circle_test ( ):

#*****************************************************************************80
#
## uniform_in_circle_test() tests uniform_in_circle().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    03 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 400

  print ( '' )
  print ( 'uniform_in_circle_test():' )
  print ( '  uniform_in_circle() returns uniform random points' )
  print ( '  inside a unit circle.' )
  print ( '' )
  print ( '  Number of points N =          ', n )

  x = uniform_in_circle ( n )

  nc = 101
  t = np.linspace ( 0.0, 2.0 * np.pi, nc )
  xc = np.zeros ( [ nc, 2 ] )
  xc[:,0] = np.cos ( t )
  xc[:,1] = np.sin ( t )

  plt.clf ( )
  plt.plot ( xc[:,0], xc[:,1], 'r', linewidth = 3 )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_circle()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_in_circle.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_ellipse ( n, A, r ):

#*****************************************************************************80
#
## uniform_in_ellipse() returns uniform points from inside an ellipse.
#
#  Discussion:
#
#    The points X in the ellipse are described by a 2 by 2
#    symmetric positive definite matrix A, and a "radius" R, such that
#
#      X' * A * X <= R^2
#
#    The upper triangular Cholesky factor of A is computed:
#
#      A = U' * U
#
#    Points Y are selected uniformly at random from the circle of radius R.
#
#      || Y || <= R
#
#    The appropriate points in the ellipse are found by solving
#
#      U * X = Y
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    04 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Reuven Rubinstein,
#    Monte Carlo Optimization, Simulation, and Sensitivity
#    of Queueing Networks,
#    Wiley, 1986.
#
#  Input:
#
#    integer N, the number of points.
#
#    real A(2,2), the matrix that describes the ellipse.
#
#    real R, the right hand side of the ellipse equation.
#
#  Output:
#
#    real P[N,2], the points.
#
  import numpy as np
#
#  Get the lower Cholesky factor L.
#
  L = np.linalg.cholesky ( A )
#
#  Get the points Y that satisfy ||Y|| <= R.
#
  p = uniform_in_circle ( n )

  p = r * p
#
#  Solve U * X = Y.
#
  for j in range ( 0, n ):
    p[j,:] = np.linalg.solve ( L.T, p[j,:] )

  return p

def uniform_in_ellipse_test ( ):

#*****************************************************************************80
#
## uniform_in_ellipse_test() tests uniform_in_ellipse().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 400
  A = np.array ( [ \
    [ 4.0, 1.0 ],
    [ 1.0, 0.5 ] ] )
  r = 1.0

  print ( '' )
  print ( 'uniform_in_ellipse_test():' )
  print ( '  uniform_in_ellipse() returns uniform random points ' )
  print ( '  from inside an ellipse.' )
  print ( '' )
  print ( '  Number of points N = ', n )
  print ( '  R                  = ', r )
  print ( '  matrix A:' )
  print ( A )

  x = uniform_in_ellipse ( n, A, r )

  plt.clf ( )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_ellipse()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_in_ellipse.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_ellipsoid ( n, A, r ):

#*****************************************************************************80
#
## uniform_in_ellipsoid() returns uniform points from inside an ellipsoid.
#
#  Discussion:
#
#    The points X in the ellipse are described by a D by D
#    symmetric positive definite matrix A, and a "radius" R, such that
#
#      X' * A * X <= R^2
#
#    The upper triangular Cholesky factor of A is computed:
#
#      A = U' * U
#
#    Points Y are selected uniformly at random from inside the sphere of radius R.
#
#      || Y || <= R
#
#    The appropriate points in the ellipse are found by solving
#
#      U * X = Y
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    04 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Reuven Rubinstein,
#    Monte Carlo Optimization, Simulation, and Sensitivity
#    of Queueing Networks,
#    Wiley, 1986.
#
#  Input:
#
#    integer N, the number of points.
#
#    real A(D,D), the matrix that describes the ellipse.
#
#    real R, the right hand side of the ellipse equation.
#
#  Output:
#
#    real P[N,2], the points.
#
  import numpy as np

  d = A.shape[0]
#
#  Get the lower Cholesky factor L.
#
  L = np.linalg.cholesky ( A )
#
#  Get the points Y that satisfy || Y || <= R.
#
  p = uniform_in_hypersphere ( n, d )

  p = r * p
#
#  Solve U * X = Y.
#
  for j in range ( 0, n ):
    p[j,:] = np.linalg.solve ( L.T, p[j,:] )

  return p

def uniform_in_ellipsoid_test ( ):

#*****************************************************************************80
#
## uniform_in_ellipsoid_test() tests uniform_in_ellipsoid().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    05 April 2022
#
#  Author:
#
#    John Burkardt
#
  from mpl_toolkits.mplot3d import Axes3D
  import matplotlib.pyplot as plt
  import numpy as np

  n = 500
  A = np.array ( [ \
    [ 4.0, 1.0, 0.0 ],
    [ 1.0, 4.0, 1.0 ],
    [ 0.0, 1.0, 0.5 ] ] )
  r = 1.0

  print ( '' )
  print ( 'uniform_in_ellipsoid_test():' )
  print ( '  uniform_in_ellipsoid() returns uniform random points ' )
  print ( '  from inside an ellipsoid x''Ax <= r^2.' )
  print ( '' )
  print ( '  Number of points N = ', n )
  print ( '  R                  = ', r )
  print ( '  matrix A:' )
  print ( A )

  p = uniform_in_ellipsoid ( n, A, r )

  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_in_ellipsoid (3d)' )
  ax1.grid ( True )
  filename = 'uniform_in_ellipsoid.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_hexagon ( n ):

#*****************************************************************************80
#
## uniform_in_hexagon() samples uniformly from the regular unit hexagon.
#
#  Discussion:
#
#    The unit hexagon has center (0,0) and "radius" 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    04 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of points.
#
#  Output:
#
#    real P(N,2), the points.
#
  from itertools import chain
  import numpy as np
#
#  Set basis vectors for each rhombus.
#
  v = np.array ( [ \
    [ -1.0,   0.0 ], \
    [  0.5,   np.sqrt ( 3.0 ) / 2.0 ], \
    [  0.5, - np.sqrt ( 3.0 ) / 2.0 ], \
    [ -1.0,   0.0 ] ] )
#
#  Assign each point randomly to rhombus 0, 1 or 2.
#
  rh = np.random.randint ( 0, 3, n )

  i = np.where ( rh == 0 )
  j = np.where ( rh == 1 )
  k = np.where ( rh == 2 )
#
#  This repulsive code is necessary because the where command returns
#  some godawful structure, not a usable list.
#
  i = tuple ( chain.from_iterable ( i ) )
  j = tuple ( chain.from_iterable ( j ) )
  k = tuple ( chain.from_iterable ( k ) )
#
#  Set barycentric coordinates of each point in its rhombus.
#
  p = np.random.rand ( n, 2 )
#
#  Convert to physical coordinates.
#
  p[i,:] = np.dot ( p[i,:], v[0:2,0:2] )
  p[j,:] = np.dot ( p[j,:], v[1:3,0:2] )
  p[k,:] = np.dot ( p[k,:], v[2:4,0:2] )

  return p

def uniform_in_hexagon_test ( ):

#*****************************************************************************80
#
## uniform_in_hexagon_test() tests uniform_in_hexagon().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 400

  print ( '' )
  print ( 'uniform_in_hexagon_test():' )
  print ( '  uniform_in_hexagon01() returns uniform random points' )
  print ( '  from inside the unit hexagon.' )
  print ( '' )
  print ( '  Number of points N =          ', n )

  x = uniform_in_hexagon ( n )

  r = np.sqrt ( 3.0 ) / 2.0

  plt.clf ( )
  plt.plot ( [1.0, 0.5, -0.5, -1.0, -0.5, 0.5, 1.0],\
             [0.0, r,    r,    0.0, -r,  -r,   0.0], 'r', linewidth = 3 )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_hexagon()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_in_hexagon.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_hypercube ( n, d ):

#*****************************************************************************80
#
## uniform_in_hypercube() creates uniform points in the unit hypercube.
#
#  Discussion:
#
#    The unit hypercube is defined as points whose components are between
#    -1 and 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    05 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of points.
#
#    integer D, the dimension of the space.
#
#  Output:
#
#    real P(N,D), the points.
#
  import numpy as np

  p = 2.0 * np.random.rand ( n, d ) - 1.0

  return p

def uniform_in_hypercube_test ( ):

#*****************************************************************************80
#
## uniform_in_hypercube_test() tests uniform_in_hypercube().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 800
  d = 3

  print ( '' )
  print ( 'uniform_in_hypercube_test():' )
  print ( '  uniform_in_hypercube() returns uniform random points ' )
  print ( '  from inside the unit hypercube in D-space.' )
  print ( '' )
  print ( '  Number of points N = ', n )
  print ( '  Space dimension D =  ', d )

  p = uniform_in_hypercube ( n, d )

  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_in_hypercube (3d)' )
  ax1.grid ( True )
  filename = 'uniform_in_hypercube.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_hypersphere ( n, d ):

#*****************************************************************************80
#
## uniform_in_hypersphere maps uniform points in the unit hypersphere.
#
#  Discussion:
#
#    The hypersphere has center 0 and radius 1.
#
#    We first generate a point ON the hypersphere, and then distribute it
#    IN the hypersphere.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    08 November 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Russell Cheng,
#    Random Variate Generation,
#    in Handbook of Simulation,
#    edited by Jerry Banks,
#    Wiley, 1998, pages 168.
#
#    Reuven Rubinstein,
#    Monte Carlo Optimization, Simulation, and Sensitivity
#    of Queueing Networks,
#    Wiley, 1986, page 232.
#
#  Input:
#
#    integer N, the number of points.
#
#    integer D, the dimension of the space.
#
#  Output:
#
#    real X(N,D), the points.
#
  import numpy as np

  exponent = 1.0 / float ( d )

  x = np.zeros ( [ n, d ] )

  for i in range ( 0, n ):
#
#  Fill a vector with normally distributed values.
#
    v = np.random.normal ( 0.0, 1.0, size = d )
#
#  Compute the length of the vector.
#
    norm = np.linalg.norm ( v )
#
#  Normalize the vector.
#
    v = v / norm
#
#  Now compute a value to map the point ON the hypersphere INTO the hypersphere.
#
    r = np.random.rand ( )

    x[i,:] = r ** exponent * v
  
  return x

def uniform_in_hypersphere_test ( ):

#*****************************************************************************80
#
## uniform_in_hypersphere_test() tests uniform_in_hypersphere().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    05 April 2022
#
#  Author:
#
#    John Burkardt
#
  from mpl_toolkits.mplot3d import Axes3D
  import matplotlib.pyplot as plt
  import numpy as np

  n = 500
  d = 3

  print ( '' )
  print ( 'uniform_in_hypersphere_test():' )
  print ( '  uniform_in_hypersphere() returns uniform random points' )
  print ( '  from inside the D-dimensional unit hypersphere.' )
  print ( '' )
  print ( '  Number of points N = ', n )
  print ( '  Space dimension D =  ', d )

  p = uniform_in_hypersphere ( n, d )

  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_in_hypersphere (3d)' )
  ax1.grid ( True )
  filename = 'uniform_in_hypersphere.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_parallelogram ( n, v1, v2, v3 ):

#*****************************************************************************80
#
## uniform_in_parallelogram() maps uniform points into a parallelogram.
#
#  Discussion:
#
#    The parallelogram is defined by three vertices, V1, V2 and V3.
#    The missing vertex V4 is equal to V2+V3-V1.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Greg Turk,
#    Generating Random Points in a Triangle,
#    in Graphics Gems,
#    edited by Andrew Glassner,
#    AP Professional, 1990, pages 24-28.
#
#  Input:
#
#    integer N, the number of points.
#
#    real V1(2), V2(2), V3(2), the vertices.
#
#  Output:
#
#    real P(N,2), the points.
#
  import numpy as np

  p = np.zeros ( [ n, 2 ] )

  for i in range ( 0, n ):

    r = np.random.rand ( 2 )

    p[i,:] = ( 1.0 - r[0] - r[1] ) * v1 \
                   + r[0]          * v2 \
                          + r[1]   * v3

  return p

def uniform_in_parallelogram_test ( ):

#*****************************************************************************80
#
## uniform_in_parallelogram_test() tests uniform_in_parallelogram().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    05 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 200
  v1 = np.array ( [ 0.75, 0.90 ] )
  v2 = np.array ( [ 0.00, 0.20 ] )
  v3 = np.array ( [ 1.50, 0.65 ] )
  v4 = v2 + v3 - v1

  print ( '' )
  print ( 'uniform_in_parallelogram_test():' )
  print ( '  uniform_in_parallelogram() returns uniform random points' )
  print ( '  from inside a parallelogram.' )
  print ( '' )
  print ( '  Number of points N =          %d' % ( n ) )
  print ( '' )
  print ( '  V1 = %12f  %12f' % ( v1[0], v1[1] ) )
  print ( '  V2 = %12f  %12f' % ( v2[0], v2[1] ) )
  print ( '  V3 = %12f  %12f' % ( v3[0], v3[1] ) )
  print ( '  V4 = %12f  %12f' % ( v4[0], v4[1] ) )

  x = uniform_in_parallelogram ( n, v1, v2, v3 )

  plt.clf ( )
  plt.plot ( [v1[0],v2[0]], [v1[1],v2[1]], 'r', linewidth = 3 )
  plt.plot ( [v2[0],v4[0]], [v2[1],v4[1]], 'r', linewidth = 3 )
  plt.plot ( [v4[0],v3[0]], [v4[1],v3[1]], 'r', linewidth = 3 )
  plt.plot ( [v3[0],v1[0]], [v3[1],v1[1]], 'r', linewidth = 3 )

  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_parallelogram()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_in_parallelogram.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_polygon ( n, nv, v ):

#*****************************************************************************80
#
## uniform_in_polygon() maps uniform points into a polygon.
#
#  Discussion:
#
#    If the polygon is regular, or convex, or at least star-shaped,
#    this routine will work.
#
#    This routine assumes that all points between the centroid and
#    any point on the boundary lie within the polygon.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of points to create.
#
#    integer NV, the number of vertices.
#
#    real V(NV,2), the vertices of the polygon, listed in
#    clockwise or counterclockwise order.
#
#  Output:
#
#    real X(N,2), the points.
#
  import numpy as np

  x = np.zeros ( [ n, 2 ] )
  t = np.zeros ( [ 3, 2 ] )
  area = np.zeros ( nv )
#
#  Find the centroid.
#
  centroid = polygon_centroid ( nv, v )
#
#  Determine the areas of each triangle.
#
  for i in range ( 0, nv ):

    if ( i < nv - 1 ):
      ip1 = i + 1
    else:
      ip1 = 0
    
    t[0,:] = v[i,:]
    t[1,:] = v[ip1,:]
    t[2,:] = centroid
    
    area[i] = triangle_area ( t )
#
#  Normalize the areas.
#
  area = area / np.sum ( area )
#
#  Replace each area by the sum of itself and all previous ones.
#
  for i in range ( 1, nv ):
    area[i] = area[i] + area[i-1]

  for j in range ( 0, n ):
#
#  Choose a triangle at random, based on areas.
#
    t = np.random.rand ( 1 )

    for k in range ( 0, nv ):

      if ( t <= area[k] ):
        i2 = k
        break
#
#  Now choose a point at random in the triangle.
#
    if ( i2 < nv - 1 ):
      i2p1 = i2 + 1
    else:
      i2p1 = 0

    r = np.random.rand ( 2 )

    if ( 1.0 < np.sum ( r ) ):
      r = 1.0 - r

    x[j,:] = ( 1.0 - r[0] - r[1] ) * v[i2,:] \
                          + r[0]   * v[i2p1,:] \
                          + r[1]   * centroid[:]

  return x

def uniform_in_polygon_test ( ):

#*****************************************************************************80
#
## uniform_in_polygon_test() tests uniform_in_polygon().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 1000
  nv = 10

  v = np.array ( [ \
    [ 0.0, 0.0 ], \
    [ 0.5, 0.3 ], \
    [ 1.0, 0.0 ], \
    [ 0.7, 0.4 ], \
    [ 1.0, 0.6 ], \
    [ 0.6, 0.6 ], \
    [ 0.5, 1.0 ], \
    [ 0.4, 0.6 ], \
    [ 0.0, 0.6 ], \
    [ 0.3, 0.4 ] ] ) 

  print ( '' )
  print ( 'uniform_in_polygon_test():' )
  print ( '  uniform_in_polygon() maps uniform' )
  print ( '  points into a polygon.' )
  print ( '' )
  print ( '  Number of points N =    ', n )
  print ( '  Number of vertices NV = ', nv )
  print ( '  Polygonal vertices:' )
  print ( v )

  p = uniform_in_polygon ( n, nv, v )
#
#  Plot.
#
  plt.clf ( )

  im1 = nv - 1
  for i in range ( 0, nv ):
    plt.plot ( [v[im1,0],v[i,0]], [v[im1,1],v[i,1]], 'r', linewidth = 3 )
    im1 = i
  plt.plot ( p[:,0], p[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_polygon()' )
  plt.grid ( True )
  plt.axis ( 'equal' )
  filename = 'uniform_in_polygon.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_sector ( n, r, t1, t2 ):

#*****************************************************************************80
#
## uniform_in_sector() maps uniform points into a circular sector.
#
#  Discussion:
#
#    The sector lies between circles with center at 0 and radius 0 and R,
#    and between rays from the center at the angles T1 and T2.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Peter Shirley,
#    Nonuniform Random Point Sets Via Warping,
#    Graphics Gems, Volume III,
#    edited by David Kirk,
#    AP Professional, 1992, 
#    ISBN: 0122861663,
#    LC: T385.G6973.
#
#  Input:
#
#    integer N, the number of points.
#
#    real R, the outer radius.
#
#    real T1, T2, the two angles, which should
#    be measured in radians, with T1 < T2.
#
#  Output:
#
#    real P(N,2), the points.
#
  import numpy as np

  u = np.random.rand ( n )
  v = np.random.rand ( n )

  t =                ( 1.0 - u ) * t1 + u * t2
  r2 = r * np.sqrt ( ( 1.0 - v ) )

  p = np.zeros ( [ n, 2 ] )
  p[:,0] = r2 * np.cos ( t )
  p[:,1] = r2 * np.sin ( t )

  return p

def uniform_in_sector_test ( ):

#*****************************************************************************80
#
## uniform_in_sector_test() tests uniform_in_sector().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 300
  r = 2.0
  t1 = 0.78
  t2 = 2.35

  print ( '' )
  print ( 'uniform_in_sector_test():' )
  print ( '  uniform_in_sector() returns uniform random points ' )
  print ( '  from inside a circular sector.' )
  print ( '' )
  print ( '  N =  ', n )
  print ( '  R  = ', r )
  print ( '  T1 = ', t1 )
  print ( '  T2 = ', t2 )

  x = uniform_in_sector ( n, r, t1, t2 )

  nc = 101
  t = np.linspace ( t1, t2, nc )

  xt = np.zeros ( [ nc, 2 ] )
  xt[:,0] = r * np.cos ( t )
  xt[:,1] = r * np.sin ( t )

  plt.clf ( )
  plt.plot ( xt[:,0], xt[:,1], 'r', linewidth = 3 )
  plt.plot ( [0.0,xt[0,0]], [0.0,xt[0,1]], 'r', linewidth = 3 )
  plt.plot ( [0.0,xt[-1,0]], [0.0,xt[-1,1]], 'r', linewidth = 3 )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_sector()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_in_sector.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_simplex ( n, d ):

#*****************************************************************************80
#
## uniform_in_simplex() maps uniform points into the unit simplex.
#
#  Discussion:
#
#    The interior of the unit D-dimensional simplex is the set of
#    points X(1:D) such that each X(I) is nonnegative, and
#    sum(X(1:M)) <= 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    07 May 2013
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Reuven Rubinstein,
#    Monte Carlo Optimization, Simulation, and Sensitivity
#    of Queueing Networks,
#    Wiley, 1986.
#
#  Input:
#
#    integer N, the number of points.
#
#    integer D, the dimension of the space.
#
#  Output:
#
#    real X(N,D), the points.
#
  import numpy as np

  p = np.zeros ( [ n, d ] )
#
#  The construction begins by sampling D+1 points from the
#  exponential distribution with parameter 1.
#
  for i in range ( 0, n ):

    e = np.random.rand ( d + 1 )

    e = - np.log ( e )

    e_sum = np.sum ( e )

    p[i,0:d] = e[0:d] / e_sum

  return p

def uniform_in_simplex_test ( ):

#*****************************************************************************80
#
## uniform_in_simplex_test() tests uniform_in_simplex().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    05 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 200
  d = 2
  v1 = np.array ( [ 0.0, 0.0 ] )
  v2 = np.array ( [ 1.0, 0.0 ] )
  v3 = np.array ( [ 0.0, 1.0 ] )

  print ( '' )
  print ( 'uniform_in_simplex_test():' )
  print ( '  uniform_in_simplex() returns uniform random points' )
  print ( '  from inside the unit simplex.' )
  print ( '' )
  print ( '  Number of points N =  ', n )
  print ( '  Spatial dimension D = ', d )

  x = uniform_in_simplex ( n, d )

  plt.clf ( )
  plt.plot ( [v1[0],v2[0]], [v1[1],v2[1]], 'r', linewidth = 3 )
  plt.plot ( [v2[0],v3[0]], [v2[1],v3[1]], 'r', linewidth = 3 )
  plt.plot ( [v3[0],v1[0]], [v3[1],v1[1]], 'r', linewidth = 3 )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_simplex()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_in_simplex.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_tetrahedron ( n, v ):

#*****************************************************************************80
#
## uniform_in_tetrahedron() returns uniform points in a tetrahedron.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Claudio Rocchini, Paolo Cignoni,
#    Generating Random Points in a Tetrahedron,
#    Journal of Graphics Tools,
#    Volume 5, Number 5, 2000, pages 9-12.
#
#  Input:
#
#    integer N, the number of points.
#
#    real V(4,3), the vertices of the tetrahedron.
#
#  Output:
#
#    real P(N,3), the points.
#
  import numpy as np

  p = np.zeros ( [ n, 3 ] );

  for i in range ( 0, n ):

    c = np.random.rand ( 3 )

    if ( 1.0 < c[0] + c[1] ):
      c[0] = 1.0 - c[0]
      c[1] = 1.0 - c[1]

    if ( 1.0 < c[1] + c[2] ):
      t = c[2]
      c[2] = 1.0 - c[0] - c[1]
      c[1] = 1.0 - t
    elif ( 1.0 < c[0] + c[1] + c[2] ):
       t = c[2]
       c[2] = c[0] + c[1] + c[2] - 1.0
       c[0] = 1.0 - c[1] - t

    c3 = 1.0 - c[0] - c[1] - c[2]

    p[i,:] = c[0] * v[0,:] \
           + c[1] * v[1,:] \
           + c[2] * v[2,:] \
           + c3   * v[3,:]

  return p

def uniform_in_tetrahedron_test ( ):

#*****************************************************************************80
#
## uniform_in_tetrahedron_test() tests uniform_in_tetrahedron().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
  from mpl_toolkits.mplot3d import Axes3D
  import matplotlib.pyplot as plt
  import numpy as np

  n = 200
  v = np.array ( [ \
    [ 1.0,  2.0,  3.0 ], \
    [ 4.0,  1.0,  2.0 ], \
    [ 2.0,  4.0,  4.0 ], \
    [ 3.0,  2.0,  5.0 ] ] )

  print ( '' )
  print ( 'uniform_in_tetrahedron_test' )
  print ( '  uniform_in_tetrahedron() returns uniform random points' )
  print ( '  from inside an arbitrary tetrahedron.' )
  print ( '' )
  print ( '  Number of points N =  ', n )
  print ( '  Vertices:' )
  print ( v )

  p = uniform_in_tetrahedron ( n, v )

  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.plot ( [v[0,0], v[1,0]], [v[0,1], v[1,1]], [v[0,2], v[1,2]], \
    'r-', linewidth = 3 )
  ax1.plot ( [v[0,0], v[2,0]], [v[0,1], v[2,1]], [v[0,2], v[2,2]], \
    'r-', linewidth = 3 )
  ax1.plot ( [v[0,0], v[3,0]], [v[0,1], v[3,1]], [v[0,2], v[3,2]], \
    'r-', linewidth = 3 )
  ax1.plot ( [v[1,0], v[2,0]], [v[1,1], v[2,1]], [v[1,2], v[2,2]], \
    'r-', linewidth = 3 )
  ax1.plot ( [v[1,0], v[3,0]], [v[1,1], v[3,1]], [v[1,2], v[3,2]], \
    'r-', linewidth = 3 )
  ax1.plot ( [v[2,0], v[3,0]], [v[2,1], v[3,1]], [v[2,2], v[3,2]], \
    'r-', linewidth = 3 )

  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_in_tetrahedron' )
  ax1.grid ( True )
  filename = 'uniform_in_tetrahedron.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_in_triangle ( n, v1, v2, v3 ):

#*****************************************************************************80
#
## uniform_in_triangle() maps uniform points into a triangle.
#
#  Discussion:
#
#    The triangle is defined by three vertices.  This routine
#    uses Turk's rule #1.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Greg Turk,
#    Generating Random Points in a Triangle,
#    in Graphics Gems,
#    edited by Andrew Glassner,
#    AP Professional, 1990, pages 24-28.
#
#  Input:
#
#    integer N, the number of points.
#
#    real V1(2), V2(2), V3(2), the vertices.
#
#  Output:
#
#    real X(N,2), the points.
#
  import numpy as np

  x = np.zeros ( [ n, 2 ] )
#
#  Generate the points using Turk's rule 1.
#
  for i in range ( 0, n ):

    r = np.random.rand ( 2 )

    a = 1.0 - np.sqrt ( r[1] )
    b = ( 1.0 - r[0] ) * np.sqrt ( r[1] )
    c = r[0] * np.sqrt ( r[1] )

    for j in range ( 0, 2 ):
      x[i,j] = ( a * v1[j] \
               + b * v2[j] \
               + c * v3[j] )

  return x

def uniform_in_triangle_test ( ):

#*****************************************************************************80
#
## uniform_in_triangle_test() tests uniform_in_triangle().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    05 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 200
  v1 = np.array ( [ 0.75, 0.90 ] )
  v2 = np.array ( [ 0.00, 0.20 ] )
  v3 = np.array ( [ 0.95, 0.65 ] )

  print ( '' )
  print ( 'uniform_in_triangle_test():' )
  print ( '  uniform_in_triangle() returns uniform random points' )
  print ( '  from inside a triangle.' )
  print ( '' )
  print ( '  Number of points N =          %d' % ( n ) )
  print ( '  V1 = %12f  %12f' % ( v1[0], v1[1] ) )
  print ( '  V2 = %12f  %12f' % ( v2[0], v2[1] ) )
  print ( '  V3 = %12f  %12f' % ( v3[0], v3[1] ) )

  x = uniform_in_triangle ( n, v1, v2, v3 )

  plt.clf ( )
  plt.plot ( [v1[0],v2[0]], [v1[1],v2[1]], 'r', linewidth = 3 )
  plt.plot ( [v2[0],v3[0]], [v2[1],v3[1]], 'r', linewidth = 3 )
  plt.plot ( [v3[0],v1[0]], [v3[1],v1[1]], 'r', linewidth = 3 )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_in_triangle()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_in_triangle.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_circle ( n ):

#*****************************************************************************80
#
## uniform_on_circle() returns randomly selected points on a circle.
#
#  Discussion:
#
#    The circle is assumed to have the equation:
#
#      x^2 + y^2 = 1
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of points to compute.
#
#  Output:
#
#    real P(N,2), points uniformly distributed on the perimeter.
#
  import numpy as np

  t = 2.0 * np.pi * np.random.rand ( n )
  p = np.zeros ( [ n, 2 ] )
  p[:,0] = np.cos ( t )
  p[:,1] = np.sin ( t )

  return p

def uniform_on_circle_test ( ):

#*****************************************************************************80
#
## uniform_on_circle_test() tests uniform_on_circle().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 100

  print ( '' )
  print ( 'uniform_on_circle_test():' )
  print ( '  uniform_on_circle() returns uniform random points ' )
  print ( '  on the surface of a circle x^2+y^2=1.' )
  print ( '' )
  print ( '  N = ', n )

  x = uniform_on_circle ( n )

  plt.clf ( )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_on_circle()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_on_circle.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_ellipse ( n, a, b ):

#*****************************************************************************80
#
## uniform_on_ellipse() returns randomly selected points on an ellipse.
#
#  Discussion:
#
#    The ellipse is assumed to have the equation:
#
#      (x/a)^2 + (y/b)^2 = 1
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    20 March 2021
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of points to compute.
#
#    real A, B: the ellipse parameters:
#
#  Output:
#
#    real P(N,2), points uniformly distributed on the perimeter of the ellipse.
#
  import numpy as np

  t = 2.0 * np.pi * np.random.rand ( n )
  p = np.zeros ( [ n, 2 ] )
  p[:,0] = a * np.cos ( t )
  p[:,1] = b * np.sin ( t )

  return p

def uniform_on_ellipse_test ( ):

#*****************************************************************************80
#
## uniform_on_ellipse_test() tests uniform_on_ellipse().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 100
  a = 2
  b = 3

  print ( '' )
  print ( 'uniform_on_ellipse_test():' )
  print ( '  uniform_on_ellipse() returns uniform random points ' )
  print ( '  on the surface of an ellipse (x/a)^2+(y/b)^2=1.' )
  print ( '' )
  print ( '  N = ', n )
  print ( '  a = ', a )
  print ( '  b = ', b )

  x = uniform_on_ellipse ( n, a, b )

  plt.clf ( )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_on_ellipse()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_on_ellipse.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_ellipsoid ( n, d, A, r ):

#*****************************************************************************80
#
## uniform_on_ellipsoid() maps uniform points onto an ellipsoid.
#
#  Discussion:
#
#    The points X on the ellipsoid are described by an N by N positive
#    definite symmetric matrix A, and a "radius" R, such that
#
#      X' * A * X = R * R
#
#    The algorithm computes the Cholesky factorization of A:
#
#      A = U' * U.
#
#    A set of uniformly random points Y is generated, satisfying:
#
#      Y' * Y = R * R.
#
#    The appropriate points in the ellipsoid are found by solving
#
#      U * X = Y
#
#    Thanks to Dr Karl-Heinz Keil for pointing out that the original
#    coding was actually correct only if A was replaced by its inverse.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    19 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Reuven Rubinstein,
#    Monte Carlo Optimization, Simulation, and Sensitivity
#    of Queueing Networks,
#    Wiley, 1986.
#
#  Input:
#
#    integer N, the number of points.
#
#    integer D, the spatial dimension.
#
#    real A(D,D), the matrix that describes the ellipsoid.
#
#    real R, the right hand side of the ellipsoid equation.
#
#  Output:
#
#    real X(N,D), the points.
#
  from scipy.linalg import cholesky
  from scipy.linalg import solve_triangular
  import numpy as np
#
#  Get the factor U such that U' * U = A.
#
  U = cholesky ( A )
#
#  Get the points Y that satisfy Y' * Y = R * R.
#
  x = uniform_on_hypersphere ( n, d )

  x = r * x
#
#  Solve U * X = Y.
#
  for i in range ( 0, n ):
    b = np.transpose ( x[i,0:d] )
    xt = solve_triangular ( U, b )
    x[i,0:d] = np.transpose ( xt )

  return x

def uniform_on_ellipsoid_test ( ):

#*****************************************************************************80
#
## uniform_on_ellipsoid_test() tests uniform_on_ellipsoid().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    19 April 2022
#
#  Author:
#
#    John Burkardt
#
  from mpl_toolkits.mplot3d import Axes3D
  import matplotlib.pyplot as plt
  import numpy as np

  n = 500
  d = 3
  A = np.array ( [ \
    [ 3.0, 1.0, 0.0 ], \
    [ 1.0, 2.0, 1.0 ], \
    [ 0.0, 1.0, 4.0 ] ] )
  r = 1.0

  print ( '' )
  print ( 'uniform_on_ellipsoid_test():' )
  print ( '  uniform_on_ellipsoid() returns uniform random points' )
  print ( '  on the surface of the D-dimensional ellipsoid x''Ax=r^2.' )
  print ( '' )
  print ( '  Number of points N = ', n )
  print ( '  Space dimension D =  ', d )
  print ( '  Radius R =           ', r )
  print ( '  Matrix A:' )
  print ( A )

  p = uniform_on_ellipsoid ( n, d, A, r )

  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_on_ellipsoid (3d)' )
  ax1.grid ( True )
  filename = 'uniform_on_ellipsoid.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_hemisphere_phong ( n, m ):

#*****************************************************************************80
#
## uniform_on_hemisphere_phong() maps uniform points onto the unit hemisphere.
#
#  Discussion:
#
#    The sphere has center 0 and radius 1.
#
#    The Phong density is used, with exponent M:
#
#    rho ( theta, phi; m ) = ( m + 1 ) * cos ( phi )^M / ( 2 * pi )
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    25 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Peter Shirley,
#    Nonuniform Random Point Sets Via Warping,
#    Graphics Gems, Volume III,
#    edited by David Kirk,
#    AP Professional, 1992, 
#    ISBN: 0122861663,
#    LC: T385.G6973.
#
#  Input:
#
#    integer N, the number of points.
#
#    integer M, the Phong exponent.
#
#  Output:
#
#    real P(N,3), the points.
#
  import numpy as np

  p = np.zeros (  [ n, 3 ] )

  power = 1.0 / float ( m + 1 )
  phi = np.random.rand ( n )
  phi = np.arccos ( ( 1.0 - phi ) ** power )

  theta = np.random.rand ( n );
  theta = 2.0 * np.pi * theta

  p[:,0] = np.cos ( theta ) * np.sin ( phi )
  p[:,1] = np.sin ( theta ) * np.sin ( phi )
  p[:,2] =                    np.cos ( phi )

  return p

def uniform_on_hemisphere_phong_test ( ):

#*****************************************************************************80
#
## uniform_on_hemisphere_phong_test() tests uniform_on_hemisphere_phong().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    25 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 500
  m = 2

  print ( '' )
  print ( 'uniform_on_hemisphere_phong_test():' )
  print ( '  uniform_on_hemisphere_phong() returns uniform random points ' )
  print ( '  on the surface of unit hemisphere with Phong density.' )
  print ( '' )
  print ( '  Number of points N = ', n )
  print ( '  Phong exponent M =   ', m )

  p = uniform_on_hemisphere_phong ( n, m )

  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_on_hemisphere_phong ' )
  ax1.grid ( True )
  filename = 'uniform_on_hemisphere_phong.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_hypercube ( n, d ):

#*****************************************************************************80
#
## uniform_on_hypercube() creates uniform points on the unit hypercube.
#
#  Discussion:
#
#    The unit hypercube is defined as points whose components are between
#    0 and 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    08 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of points.
#
#    integer D, the dimension of the space.
#
#  Output:
#
#    real P(N,D), the points.
#
  import numpy as np
#
#  Choose random points within the cube of radius 1.
#
  p = 2.0 * np.random.rand ( n, d ) - 1.0
#
#  For each point, select a coordinate at random, and set it to +1 or -1.
#
  for i in range ( 0, n ):

    j = np.random.random_integers ( 0, d - 1, 1 )
    s = 2.0 * np.random.random_integers ( 0, 1, 1 ) - 1.0
    p[i,j] = s

  return p

def uniform_on_hypercube_test ( ):

#*****************************************************************************80
#
## uniform_on_hypercube_test() tests uniform_on_hypercube().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    08 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 800
  d = 3

  print ( '' )
  print ( 'uniform_on_hypercube_test():' )
  print ( '  uniform_on_hypercube() returns uniform random points ' )
  print ( '  on the surface of the unit hypercube in D-space.' )
  print ( '' )
  print ( '  Number of points N = ', n )
  print ( '  Space dimension D =  ', d )

  p = uniform_on_hypercube ( n, d )

  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_on_hypercube (3d)' )
  ax1.grid ( True )
  filename = 'uniform_on_hypercube.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_hypersphere ( n, d ):

#*****************************************************************************80
#
## uniform_on_hypersphere() samples surface of the unit hypersphere.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    03 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Russell Cheng,
#    Random Variate Generation,
#    in Handbook of Simulation,
#    edited by Jerry Banks,
#    Wiley, 1998, pages 168.
#
#    George Marsaglia,
#    Choosing a point from the surface of a sphere,
#    Annals of Mathematical Statistics,
#    Volume 43, Number 2, April 1972, pages 645-646.
#
#    Reuven Rubinstein,
#    Monte Carlo Optimization, Simulation, and Sensitivity
#    of Queueing Networks,
#    Wiley, 1986, page 234.
#
#  Input:
#
#    integer N, the number of points.
#
#    integer D, the dimension of the space.
#
#  Output:
#
#    real X(N,D), the points.
#
  import numpy as np

  x = np.random.normal ( 0.0, 1.0, size = ( n, d ) )
  v = np.linalg.norm ( x, axis = 1 )
  for j in range ( 0, d ):
    x[:,j] = x[:,j] / v

  return x

def uniform_on_hypersphere_test ( ):

#*****************************************************************************80
#
## uniform_on_hypersphere_test() tests uniform_on_hypersphere().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    06 April 2022
#
#  Author:
#
#    John Burkardt
#
  from mpl_toolkits.mplot3d import Axes3D
  import matplotlib.pyplot as plt
  import numpy as np

  n = 200
  d = 3

  print ( '' )
  print ( 'uniform_on_hypersphere_test():' )
  print ( '  uniform_on_hypersphere() computes points uniformly distributed' )
  print ( '  on the surface of the unit hypersphere.' )
  print ( '' )
  print ( '  Number of points N = ', n )
  print ( '  Space dimension D =  ', d )

  p = uniform_on_hypersphere ( n, d )

  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_on_hypersphere (3d)' )
  ax1.grid ( True )
  filename = 'uniform_on_hypersphere.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_simplex ( n, d ):

#*****************************************************************************80
#
## uniform_on_simplex() maps uniform points onto the unit simplex.
#
#  Discussion:
#
#    The surface of the unit D-dimensional simplex is the set of points
#    X(1:D) such that each X(I) is nonnegative,
#    every X(I) is no greater than 1, and
#
#    ( X(I) = 0 for some I, or sum ( X(1:D) ) = 1. )
#
#    In D dimensions, there are D sides, and one main face.
#    This code picks a point uniformly with respect to "area".
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    08 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Reuven Rubinstein,
#    Monte Carlo Optimization, Simulation, and Sensitivity
#    of Queueing Networks,
#    Wiley, 1986.
#
#  Input:
#
#    integer N, the number of points.
#
#    integer D, the dimension of the space.
#
#  Output:
#
#    real X(N,D), the points.
#
  import numpy as np

  p = np.zeros ( [ n, d ] );
#
#  The construction begins by sampling M points from the
#  exponential distribution with parameter 1.
#
  for i in range ( 0, n ):

    e = np.random.rand ( d )
    e = - np.log ( e )
    e_sum = np.sum ( e )

    p[i,:] = e / e_sum
#
#  The point X is now on the "main" face of the unit simplex.
#
#  Based on their relative areas, choose a side of the simplex,
#  or the main face.
#
    area1 = np.sqrt ( d )
    area2 = d

    r = np.random.rand ( 1 )
#
#  If we choose to move the point from the main face,
#  set a random coordinate to 0.
#
    if ( area1 / ( area1 + area2 ) < r ):
      j = np.random.randint ( 0, d, 1 )
      p[i,j] = 0.0

  return p

def uniform_on_simplex_test ( ):

#*****************************************************************************80
#
## uniform_on_simplex_test() tests uniform_on_simplex().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    08 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 200
  d = 2
  v1 = np.array ( [ 0.0, 0.0 ] )
  v2 = np.array ( [ 1.0, 0.0 ] )
  v3 = np.array ( [ 0.0, 1.0 ] )

  print ( '' )
  print ( 'uniform_on_simplex_test():' )
  print ( '  uniform_on_simplex() returns uniform random points' )
  print ( '  from inside the unit simplex.' )
  print ( '' )
  print ( '  Number of points N =  ', n )
  print ( '  Spatial dimension D = ', d )

  x = uniform_on_simplex ( n, d )

  plt.clf ( )
  plt.plot ( [v1[0],v2[0]], [v1[1],v2[1]], 'r', linewidth = 3 )
  plt.plot ( [v2[0],v3[0]], [v2[1],v3[1]], 'r', linewidth = 3 )
  plt.plot ( [v3[0],v1[0]], [v3[1],v1[1]], 'r', linewidth = 3 )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_on_simplex()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_on_simplex.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_sphere_patch ( n, phi1, phi2, theta1, theta2 ):

#*****************************************************************************80
#
## uniform_on_sphere_patch() maps uniform points to a spherical patch.
#
#  Discussion:
#
#    The sphere has center 0 and radius 1.
#
#    A spherical patch on the surface of the unit sphere contains those 
#    points with radius R = 1 and angles (THETA,PHI) such that
#
#      0.0 <= THETA1 <= THETA <= THETA2 <= 2 * PI
#      0.0 <= PHI1   <= PHI   <= PHI2   <=     PI
#
#    mapped to:
#
#      x = cos ( theta ) * sin ( phi )
#      y = sin ( theta ) * sin ( phi )
#      z =                 cos ( phi )
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    20 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Peter Shirley,
#    Nonuniform Random Point Sets Via Warping,
#    Graphics Gems, Volume III,
#    edited by David Kirk,
#    AP Professional, 1992, 
#    ISBN: 0122861663,
#    LC: T385.G6973.
#
#  Input:
#
#    integer N, the number of points.
#
#    real PHI1, PHI2, the latitudinal angle range.
#
#    real THETA1, THETA2, the longitudinal angle range.
#
#  Output:
#
#    real P(N,3), the points, in (X,Y,Z) coordinates.
#
  import numpy as np

  r = np.random.rand ( n )
  theta = ( 1.0 - r ) * theta1 + r * theta2

  r = np.random.rand ( n )
  phi = np.arccos ( ( 1.0 - r ) * np.cos ( phi1 ) + r * np.cos ( phi2 ) )

  p = np.zeros ( [ n, 3 ] )
  p[:,0] = np.cos ( theta ) * np.sin ( phi )
  p[:,1] = np.sin ( theta ) * np.sin ( phi )
  p[:,2] =                    np.cos ( phi )

  return p

def uniform_on_sphere_patch_test ( ):

#*****************************************************************************80
#
## uniform_on_sphere_patch_test() tests uniform_on_sphere_patch().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    20 April 2022
#
#  Author:
#
#    John Burkardt
#
  from mpl_toolkits.mplot3d import Axes3D
  import matplotlib.pyplot as plt
  import numpy as np

  n = 500
  phi1 =    80.0 * ( np.pi / 180.0 )
  phi2 =   100.0 * ( np.pi / 180.0 )
  theta1 =   0.0 * ( np.pi / 180.0 )
  theta2 = 360.0 * ( np.pi / 180.0 )

  print ( '' )
  print ( 'uniform_on_sphere_patch_test():' )
  print ( '  uniform_on_sphere_patch() maps uniform' )
  print ( '  points onto a patch of the unit sphere.' )
  print ( '  The points are returned as (X,Y,Z) coordinates.' )
  print ( '' )
  print ( '  Number of points N =          ', n )
  print ( '  Latitudinal angle PHI1 =      ', phi1 )
  print ( '  Latitudinal angle PHI2 =      ', phi2 )
  print ( '  Longitudinal angle THETA1 =   ', theta1 )
  print ( '  Longitudinal angle THETA2 =   ', theta2 )

  p = uniform_on_sphere_patch ( n, phi1, phi2, theta1, theta2 )
#
#  Make a plot.
#
  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_on_sphere_patch()' )
  ax1.grid ( True )
  filename = 'uniform_on_sphere_patch.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_sphere_triangle ( n, v1, v2, v3 ):

#*****************************************************************************80
#
## uniform_on_sphere_triangle(): sample spherical triangle.
#
#  Discussion:
#
#    The sphere has center 0 and radius 1.
#
#    A spherical triangle on the surface of the unit sphere contains those 
#    points with radius R = 1, bounded by the vertices V1, V2, V3.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    24 August 2010
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    James Arvo,
#    Stratified sampling of spherical triangles,
#    Computer Graphics Proceedings, Annual Conference Series, 
#    ACM SIGGRAPH '95, pages 437-438, 1995.
#
#  Input:
#
#    integer N, the number of points.
#
#    real V1(1,3), V2(1,3), V3(1,3), the XYZ coordinates of
#    the vertices of the spherical triangle.
#
#  Output:
#
#    real P(N,3), the XYZ coordinates of the sample points.
#
  import numpy as np

  p = np.zeros ( [ n, 3 ] )
#
#  Compute the sides, angles, and area of the spherical triangle
#  for now, we assume R = 1.
#
  r = 1.0

  a, b, c = sphere_triangle_vertices_to_sides ( r, v1, v2, v3 )

  alpha, beta, gamma = sphere_triangle_sides_to_angles ( r, a, b, c )

  area = sphere_triangle_angles_to_area ( r, alpha, beta, gamma )

  for i in range ( 0, n ):
#
#  Select the new area.
#
    xsi1 = np.random.rand ( 1 )
    area_hat = xsi1 * area
#
#  Compute the sine and cosine of the angle phi.
#
    s = np.sin ( area_hat - alpha )
    t = np.cos ( area_hat - alpha )
#
#  Compute the pair that determines beta_hat.
#
    u = t - np.cos ( alpha )
    v = s + np.sin ( alpha ) * np.cos ( c )
#
#  Q is the cosine of the new edge length b_hat.
#
    q = ( ( v * t - u * s ) * np.cos ( alpha ) - v ) \
      / ( ( v * s + u * t ) * np.sin ( alpha ) )
#
#  V31 = normalized ( V3 - ( V3 dot V1 ) * V1 )
#
    w = np.dot ( v3, v1 )
    v31 = v3 - w * v1
    v31 = v31 / np.linalg.norm ( v31 )
#
#  V4 is the third vertex of the subtriangle V1, V2, V4.
#
    v4 = q * v1 + np.sqrt ( 1.0 - q * q ) * v31
#
#  Select cos theta, which will sample along the edge from V2 to V4.
#
    xsi2 = np.random.rand ( 1 )
    z = 1.0 - xsi2 * ( 1.0 - np.dot ( v4, v2 ) )
#
#  V42 = normalized ( V4 - ( V4 dot V2 ) * V2 )
#
    w = v4 * v2
    v42 = v4 - w * v2
    v42 = v42 / np.linalg.norm ( v42 )
#
#  Construct the point.
#
    p[i,:] = z * v2 + np.sqrt ( 1.0 - z * z ) * v42

  return p

def uniform_on_sphere_triangle_test ( ):

#*****************************************************************************80
#
## uniform_on_sphere_triangle_test() tests uniform_on_sphere_triangle().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    25 April 2022
#
#  Author:
#
#    John Burkardt
#
  from mpl_toolkits.mplot3d import Axes3D
  import matplotlib.pyplot as plt
  import numpy as np

  n = 500
  v1 = np.array ( [ 1.0, 0.0, 0.0 ] )
  v2 = np.array ( [ 0.0, 1.0, 0.0 ] )
  v3 = np.array ( [ 0.0, 0.0, 1.0 ] )

  print ( '' )
  print ( 'uniform_on_sphere_triangle_test():' )
  print ( '  uniform_on_sphere_triangle() maps uniform' )
  print ( '  points onto a triangle of the unit sphere.' )
  print ( '  The points are returned as (X,Y,Z) coordinates.' )
  print ( '' )
  print ( '  Number of points N =          ', n )
  print ( '  Vertex V1 =      ', v1 )
  print ( '  Vertex V2 =      ', v2 )
  print ( '  Vertex V3 =      ', v3 )

  p = uniform_on_sphere_triangle ( n, v1, v2, v3 )
#
#  Make a plot.
#
  fig1 = plt.figure ( )
  ax1 = fig1.add_subplot ( 111, projection = '3d' )
  ax1.scatter ( p[:,0], p[:,1], p[:,2], 'bo' );
  ax1.set_xlabel ( '--X axis--' )
  ax1.set_ylabel ( '--Y axis--' )
  ax1.set_zlabel ( '--Z axis--' )
  ax1.set_title ( 'uniform_on_sphere_triangle()' )
  ax1.grid ( True )
  filename = 'uniform_on_sphere_triangle.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_on_triangle ( n, v ):

#*****************************************************************************80
#
## uniform_on_triangle() maps uniform points onto the boundary of a triangle.
#
#  Discussion:
#
#    The triangle is defined by the three vertices V1, V2, V3.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of points.
#
#    real V(3,2), the vertices of the triangle.
#
#  Output:
#
#    real P(N,2), the points.
#
  import numpy as np

  p = np.zeros ( [ n, 2 ] );

  l1 = np.sqrt ( ( v[1,0] - v[0,0] ) ** 2 \
               + ( v[1,1] - v[0,1] ) ** 2 )

  l2 = np.sqrt ( ( v[2,0] - v[1,0] ) ** 2 \
               + ( v[2,1] - v[1,1] ) ** 2 )

  l3 = np.sqrt ( ( v[0,0] - v[2,0] ) ** 2 \
               + ( v[0,1] - v[2,1] ) ** 2 )
  
  for i in range ( 0, n ):
#
#  R can be regarded as the distance of the point on the perimeter,
#  as measured from the origin, along the perimeter.
#
    r = np.random.rand ( 1 )

    r = ( l1 + l2 + l3 ) * r
#
#  Case 1: between V1 and V2.
#
    if ( r <= l1 ):

      s = ( l1 - r ) / l1
      t =        r   / l1
      p[i,0] = s * v[0,0] + t * v[1,0]
      p[i,1] = s * v[0,1] + t * v[1,1]
#
#  Case 2: between V2 and V3.
#
    elif ( r <= l1 + l2 ):

      s = ( l2 - r + l1 ) / l2
      t = (      r - l1 ) / l2
      p[i,0] = s * v[1,0] + t * v[2,0]
      p[i,1] = s * v[1,1] + t * v[2,1]
#
#  Case 3: between V3 and V1.
#
    else:

      s = ( l3 - r + l1 + l2 ) / l3
      t = (      r - l1 - l2 ) / l3
      p[i,0] = s * v[2,0] + t * v[0,0]
      p[i,1] = s * v[2,1] + t * v[0,1]

  return p

def uniform_on_triangle_test ( ):

#*****************************************************************************80
#
## uniform_on_triangle_test() tests uniform_on_triangle().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    05 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 100
  v = np.array ( [ \
    [ 0.75, 0.90 ], \
    [ 0.00, 0.20 ], \
    [ 0.95, 0.65 ] ] )

  print ( '' )
  print ( 'uniform_on_triangle_test():' )
  print ( '  uniform_in_triangle() returns uniform random points' )
  print ( '  on the surface of a triangle.' )
  print ( '' )
  print ( '  Number of points N =          %d' % ( n ) )
  print ( '  Triangle vertices:' )
  print ( v )

  x = uniform_on_triangle ( n, v )
#
#  Plot.
#
  plt.clf ( )
  plt.plot ( [v[0,0],v[1,0]], [v[0,1],v[1,1]], 'r', linewidth = 3 )
  plt.plot ( [v[1,0],v[2,0]], [v[1,1],v[2,1]], 'r', linewidth = 3 )
  plt.plot ( [v[2,0],v[0,0]], [v[2,1],v[0,1]], 'r', linewidth = 3 )
  plt.plot ( x[:,0], x[:,1], 'bo' )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_on_triangle()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_on_triangle.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def uniform_walk ( n, d ):

#*****************************************************************************80
#
## uniform_walk() generates points on a uniform random walk.
#
#  Discussion:
#
#    The first point is at the origin.  Uniform random numbers are
#    generated to determine the direction of the next step, which
#    is always of length 1, and in coordinate direction.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    22 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of points.
#
#    integer D, the dimension of the space.
#
#  Output:
#
#    real P(N,D), the points.
#
  import numpy as np

  p = np.zeros ( [ n, d ] )
#
#  Choose random directions of +1 and -1.
#
  dir = np.random.randint ( 0, high = d, size = n - 1 )
  sgn = np.random.rand ( n - 1 ) - 0.5

  for i in range ( 1, n ):
#
#  Start at the previous point.
#
    p[i,:] = p[i-1,:]
#
#  Pick which direction to alter.
#
    j = dir[i-1]
#
#  Move +1 or -1 in that direction.
#
    if ( sgn[i-1] < 0.0 ):
      p[i,j] = p[i,j] - 1.0
    else:
      p[i,j] = p[i,j] + 1.0

  return p

def uniform_walk_test ( ):

#*****************************************************************************80
#
## uniform_walk_test() tests uniform_walk().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    22 April 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  n = 100
  d = 2

  print ( '' )
  print ( 'uniform_walk_test():' )
  print ( '  uniform_walk() generates points on a uniform random walk' )
  print ( '  of fixed steps in coordinate directions.' )
  print ( '' )
  print ( '  Number of points N =  %d' % ( n ) )
  print ( '  Spatial dimension D = %d' % ( d ) )

  p = uniform_walk ( n, d )
#
#  Plot.
#
  plt.clf ( )
  plt.plot ( p[:,0], p[:,1], 'k-' )
  plt.plot ( p[:,0], p[:,1], 'b.' )
  plt.plot ( p[0,0], p[0,1], 'g.', markersize = 15 )
  plt.plot ( p[n-1,0], p[n-1,1], 'r.', markersize = 15 )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'uniform_walk()' )
  plt.grid ( 'True' )
  plt.axis ( 'equal' )
  filename = 'uniform_walk.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  random_data_test ( )
  timestamp ( )