#! /usr/bin/env python3
#
def motion_3d ( ):

#*****************************************************************************80
#
## brownian_motion_simulation_test() tests brownian_motion_simulation().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    22 February 2025
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'motion_3d():' )
  print ( '  Simulate brownian_motion_in 3D.' )

  m = 3
  n = 1001
  d = 10.0
  t = 1.0
  x = brownian_motion_simulation ( m, n, d, t )
  motion_3d_display ( m, n, x )

  return

def brownian_motion_simulation ( m, n, d, t ):

#*****************************************************************************80
#
## brownian_motion_simulation() simulates Brownian motion.
#
#  Discussion:
#
#    Thanks to Feifei Xu for pointing out a missing factor of 2 in the
#    stepsize calculation, 08 March 2016.
#
#    Thanks to Joerg Peter Pfannmoeller for pointing out a missing factor
#    of M in the stepsize calculation, 23 April 2018.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    02 January 2024
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer M, the spatial dimension.
#
#    integer N, the number of time steps to take, plus 1. 
#
#    real D, the diffusion coefficient.  
#
#    real T, the total time.
#
#  Output:
#
#    real X(N,M), N successive locations of a particle in M dimensional space.
#
  import numpy as np
#
#  Set the time step.
#
  dt = t / float ( n - 1 )
#
#  Choose the standard deviation.
#
  sigma = np.sqrt ( 2.0 * m * d * dt )
#
#  Make space for the data.
#
  x = np.zeros ( [ n, m ] )
#
#  Compute the individual steps.
#
  for j in range ( 1, n ):
#
#  S is the stepsize
#
    s = sigma * np.random.standard_normal ( )
#
#  Direction is random.
#
    u = np.random.standard_normal ( m )
    u = u / np.linalg.norm ( u )
#
#  Each position is the sum of the previous steps.
#
    x[j,0:m] = x[j-1,0:m] + s * u[0:m]

  return x

def motion_3d_display ( m, n, x ):

#*****************************************************************************80
#
## brownian_motion_display() displays successive Brownian motion positions.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    26 April 2018
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer M, the spatial dimension.
#    M should be 1, 2 or 3.
#
#    integer N, the number of time steps. 
#
#    real X(N,M), the particle positions.
#
  import matplotlib.pyplot as plt
  import numpy as np
  from mpl_toolkits.mplot3d import Axes3D

  fig = plt.figure ( )
  ax = fig.add_subplot ( 111, projection = '3d' )
  ax.plot ( x[:,0], x[:,1], x[:,2], 'b', linewidth = 2 )
  ax.scatter ( x[0,0], x[0,1], x[0,2], c = 'g', marker = 'o', s = 100 )
  ax.scatter ( x[n-1,0], x[n-1,1], x[n-1,2], c = 'r', marker = 'o', s = 100 ) 
  ax.grid ( True )
  ax.set_xlabel ( '<--X-->' )
  ax.set_ylabel ( '<--Y-->' )
  ax.set_zlabel ( '<--Z-->' )
  plt.title ( 'Brownian motion simulation in 3D' )
  filename = 'motion_3d.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )

  return

if ( __name__ == '__main__' ):
  motion_3d()