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

#*****************************************************************************80
#
## motion_1d() simulates 1D Brownian motion
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    22 February 2025
#
#  Author:
#
#    John Burkardt
#
  import matplotlib
  import numpy as np

  print ( '' )
  print ( 'motion_1d():' )
  print ( '  Simulation 1D Brownian motion.' )

  m = 1
  n = 101
  d = 10.0
  t = 1.0
  x = brownian_motion_simulation ( m, n, d, t )
  motion_1d_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 )
#
#  Make space for the data.
#
  x = np.zeros ( [ n, m ] )
#
#  Define the standard deviation.
#
  sigma = np.sqrt ( 2.0 * m * d * dt )
#
#  Compute the individual steps.
#
  for j in range ( 1, n ):
#
#  S is the stepsize
#
    s = np.random.standard_normal ( )
#
#  Direction is unit random vector.
#
    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_1d_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

  y = np.linspace ( 0, n - 1, n ) / float ( n - 1 )
  plt.plot ( x[:,0], y[:], 'b', linewidth = 2 )
  plt.plot ( x[0,0], y[0], 'g.', markersize = 35 )
  plt.plot ( x[n-1,0], y[n-1], 'r.', markersize = 35 )
  plt.grid ( True )
  plt.xlabel ( '<--X-->' )
  plt.ylabel ( '<--Time-->' )
  plt.title ( 'Brownian motion simulation in 1D' )
  filename = 'motion_1d.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )

  return

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