#! /usr/bin/env python3
#
def two_flies_circle ( n ):

#*****************************************************************************80
#
## two_flies_circle(): distance between two random flies on unit circle.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    09 July 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer n: the number of observations.
#
#  Output:
#
#    real d[n]: the observed distance values.
#
  from numpy.random import default_rng
  import numpy as np

  rng = default_rng ( )

  d = np.zeros ( n )

  for i in range ( 0, n ):

    [ r1, r2 ] = rng.random ( 2 )
    r = np.sqrt ( r1 )
    theta = 2.0 * np.pi * r2
    x1 = r * np.cos ( theta )
    y1 = r * np.sin ( theta )

    [ r1, r2 ] = rng.random ( 2 )
    r = np.sqrt ( r1 )
    theta = 2.0 * np.pi * r2
    x2 = r * np.cos ( theta )
    y2 = r * np.sin ( theta )

    d[i] = np.sqrt ( ( x1 - x2 )**2 + ( y1 - y2 )**2 )

  return d

def two_flies_circle_test ( ):

#*****************************************************************************80
#
## two_flies_circle_test() tests two_flies_circle().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    08 July 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  mean_exact = 128.0 / 45.0 / np.pi

  print ( '' )
  print ( 'two_flies_circle_test()' )
  print ( '  Two flies land at random on a circular plate of radius 1.' )
  print ( '  What is the average distance between them?' )
  print ( '  Exact average is', mean_exact )
  print ( '' )
  print ( '      n    dmin   dmean    dmax    dvar    dstd  ||dmean-exact||' )
  print ( '' )

  for n in [ 10, 100, 1000, 10000 ]:

    d = two_flies_circle ( n )

    d_min = np.min ( d )
    d_mean = np.mean ( d )
    d_max = np.max ( d )
    d_var = np.var ( d )
    d_std = np.std ( d )

    print ( '  %5d  %6.4f  %6.4f  %6.4f  %6.4f  %6.4f  %6.4f' \
    % ( n, d_min, d_mean, d_max, d_var, d_std, np.abs ( d_mean - mean_exact ) ) )

  return

def timestamp ( ):

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

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

  return

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